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SECTION  I 


GENERAL  DESCRIPTION 

A.  Overview 

The  purpose  of  this  report  is  to  document  a  series  of  four  computer 
programs  that  are  used  to  compute  the  flowfield  within  a  two-dimensional 
mixed-compression  high  speed  aircraft  inlet.  In  Section  I,  a  brief 
description  is  presented  of  the  physical  problem  and  the  mathematical 
model.  The  numerical  methods  are  discussed  in  Section  II,  and  the 
limitations  of  the  approach  are  indicated.  A  brief  summary  of  the  four 
programs  is  presented,  together  with  the  general  sequence  of  application. 

In  Section  III,  the  coordinate  system  programs  are  discussed  in  detail, 
with  emphasis  on  the  pertinent  criteria  for  successful  implementation. 

In  Section  IV,  the  details  of  the  Navier-Stokes  code  employed  for 
solution  of  the  inlet  flowfield  are  presented.  Tn  addition,  the  details 
of  a  simple  utility  program  used  to  interpolate  flowfield  data  are 
discussed.  In  each  section,  examples  of  calculations  performed  on  the 
CYBER  74  and  175  computers  at  WPAFB  are  presented  to  illustrate  application 
of  the  numerical  codes. 

B.  Description  of  Physical  Problem 

The  basic  problem  is  the  computation  of  the  steady  flowfield  within  a 
two-dimensional  high  speed  inlet.  As  indicated  in  Figure  1,  the  incoming 
supersonic  flow  is  deflected  by  a  pattern  of  oblique  shock  waves  formed  bv 
the  general  curvilinear  shape  of  the  ramp  and  cowl  surfaces.  The  boundary 
layers  are  turbulent  over  nearly  the  entire  length  of  the  inlet,  owing  to 
the  typically  high  Reynolds  numbers.  Boundary  layer  bleed  is  distributed 
along  the  walls  in  order  to  prevent  flew  separation  at  the  interaction  of 


the  shock  waves  and  the  boundary  layers  on  the  ramp  and  cowl.  A  terminal 
shock  may  be  positioned  near  the  inlet  throat. 

The  computer  programs  described  herein  determine  the  flowfield  within 
the  inlet,  including  both  inviscid  and  viscous  regions,  by  integration 
of  the  full  mean  compressible  Navier-Stokes  equations.  Turbulence  is 
represented  by  an  algebraic  turbulent  eddy  viscosity. 

C.  Description  of  Mathematical  Model 

1 .  Coordinate  Transformation 

In  order  to  handle  the  general  curvilinear  inlet  geometry,  a  numerical 
coordinate  transformation  is  employed.  The  purpose  of  the  transformation 
is  to  provide  a  set  of  curvilinear  coordinates  ?(x,y)  and  n(x,y)  that 
are  contoured  to  the  inlet  shape.  As  indicated  schematically  in  Figure  2, 
the  cow],  and  ramp  surfaces  are  taken  to  be  coincident  with  portions  of 
the  contours  n(x,y)  =  0  and  q(x,y)  =  1,  respectively.  The  upstream  and 
downstream  boundaries  are  defined  by  the  lines  c(x,y)  =  0  and  t;(x,y)  =  1, 
respectively.  The  n  coordinate  increases  in  a  general  cross-stream 
direction  across  the  inlet,  while  the  ^  coordinate  increases  basically  in 
the  streamwise  direction.  The  clear  advantages  of  a  general  coordinate 
transformation  of  this  type  are  the  simplicity  of  application  of  the  fluid 
dynamical  boundary  conditions  (e.g.,  adiabatic  wall,  no  slip,  mass  bleed) 
for  the  ramp  and  cowl  and  the  inherent  accuracy  of  a  surface-oriented 
coordinate  system  in  the  boundary  layers.  By  construction,  the  coordinate 
transformation  maps  the  desired  domain  of  the  inlet  flowfield  in  the 
physical  or  x-y  plane  into  the  unit  square  in  the  transformed  or  C-n 
plane.  For  example,  the  portion  of  the  n  =  0  line  that  corresponds  to 
the  cowl  surface  is  mapped  into  a  segment  of  the  lower  boundary  of  the 


unit  square  in  the  transformed  plane,  where  application  ol  the  appropriate 
boundary  conditions  is  facilitated. 

The  coordinates  5(x,y)  and  n(x,y)  are  obtained  using  the  basic 
approach  of  Thompson. ^  They  are  taken  to  satisfy  the  following  equations: 


2 

r?  =  o 


2  2  2  2  2 

where  V  is  the  Laplacean  operator  9  /fx’’  +  T~/3y  .  The  coordinates 

C(x,v)  and  ri(x,y)  are  subject  to  Dirichlet  houndarv  conditions  as 

illustrated  in  Figure  2  (e.g.,  5  =  0  on  the  upstream  boundary) .  The 

coordinate  transformation  is  generally  non-orthcgonal  .  The  source  terms 

on  the  right-hand  side  of  (2)  are  employed  to  control  the  nosh  spacing 

in  the  n-direction  in  order  to  provide  accurate  resolution  within  the 

houndarv  layers  on  the  ramp  and  cowl.  The  expression  employed  for 

T(^,n)  is 


T(c,n) 


'Cl/nl 


H-C 0 / ( 1  -  n2) 


(2) 


where  C^,  C9,  and  p 9  are  slowly  varying  functions  o!  the  streanwise 
variable  r,  .  As  discussed  later  in  Section  Ill,  these  qn  nit  i  t  i  es  are 
determined  by  the  requirements  of  accurate  resol ution  of  the  boundary 
layers  and  controllable  mesh  spacing  near  the  wuiis.  The  patticularlv 
simple  form  of  the  source  term  in  (2)  provides  for  subs! ant i a l  automation 
of  the  coordinate  generation  process.  A  different  tormul  1 1 ion  for  the 
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source  terms  has  been  employed  by  Thompson  utilizing  a  series  of 
exponential  functions.  Although  possessing  greater  generality  than  the 
present  technique,  its  operation  typically  requires  a  trial-and-error 
procedure  in  order  to  generate  a  suitable  coordinate  system.  In  contrast, 
the  present  technique  permits  direct  generation  of  a  suitable  coordinate 
transformation  satisfying  the  necessary  requirements  of  accurate  resolution 
and  controllable  mesh  spacing  within  the  boundary  layers. 

The  governing  fluid  dynamic  equations  are  solved  in  the  transformed 
?-n  plane.  From  knowledge  of  the  inverse  transformation  x(t,n)  and 
y(?,n),  the  flowfield  within  the  inlet  in  the  physical  plane  is  obtained. 

2 .  Navier-Stokes  Equations 

The  governing  equations  are  the  full  mean  compressible  Navier-Stokes 

2 

equations  utilizing  mass-averaged  variables  for  two-dimensional  turbulent 

flow.  Written  in  strong  conservation  form  utilizing  the  transformed 
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coordinates  £(x,y)  and  n(x,y),  the  equations  are  ' 


where 
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at  a c  an 


'  p ' 


pu 

pv 

pe 


F-j< 


PU 


puU  +  C  ( p  -  T  )  -  c  T 

x  xx  y  xy 

pvU  +  £  (p  -  t  )-CT 

y  K  yy  x  xy 

(pe  +  p)U  +  r,  p  +  t  p, 
xx  y  y 
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e  is  the  turbulent  eddy  viscosity,  XT  =  -  (2/3) (p  +  e)  and  div  v  =  3u/ 3x  +  3v/ 3y . 
The  partial  derivatives  in  (8)  are  replaced  by  derivatives  with  respect  to 
the  transformed  variables  £  and  n  by  means  of  the  chain  rule  (e.g., 

3u/3x  =  ?x(3u/30  +  nx(3u/3n)).  The  quantities  (?x  and  Ry  in  (5)  and  (6) 
are 
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where  Q^,  are  components  of  the  heat  flux  given  by 


Q  e 


_T_j!fi 

Pr  Pr  J  3x 
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-y 


y  +  e 


Pr  PrJ  3y 


3e . 
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(10) 


where  y  =  cp/cv  is  the  ratio  of  specific  heats  (taken  to  be  1.4)  and  Pr 
and  Prt  are  the  molecular  and  turbulent  Prandtl  numbers,  respectively, 
with  values  of  0.72  and  0.90.  Although  Pr^  is,  in  general,  a  function 
of  position  within  the  boundary  layer,  numerical  computations  of  high 
speed  turbulent  boundary  layers  have  shown  only  a  weak  dependence  on 
Pr  variations.^ 

The  turbulent  eddy  viscosity  z  is  given  by  the  tv.o-layer  equilibrium 
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eddy  viscosity  model  of  Cebeci-Smith  (except  as  noted  later)  with  the 

9 

transition  model  of  Dhawan  and  Narasimha.  Within  the  inner  region,  the 

eq^  is  given  by 


equilibrium  eddy  viscosity  z 


(12) 


cq. 


pk2Uref  6i  F(s) 


where  n  =  distance  normal  to  surface 

s  =  distance  along  surface  from  leading  edge 
k  =  0.40  (von  Karman's  constant) 
k2  =  0.0168 

r(s)  =  transition  factor 

U  ^  =  mean  velocity  outside  boundary  layer 

U  =  component  of  velocity  parallel  to  the  surface 


c 


(1  -  u/U  f)  dn 


6  =  local  boundary  layer  thickness 
and  D  is  the  modified  Van  Driest  damping  factor  given  by 


n  /  t  p 

D  =  1_exp( - ^-N) 


where  N 


a  modification  due  to  mass  bleed: 


(13) 


N 


exp 


5.9 


p  m 
w 


(14) 


In  the  above,  the  subscript  w  implies  evaluation  at  the  wall,  with  t 

w 

denoting  the  wall  shear  stress,  and  m  is  the  normal  mass  flux  at  the 

■  a  A 

surface  (i.e.,  m  =  pv  •  n,  where  n  is  the  outward  normal  at  the  wall, 

pointing  into  the  fluid) .  The  above  expression  for  the  Van  Driest 

£ 

damping  factor  D  differs  from  that  of  Cebeci-Smith  in  two  respects. 
First,  the  density  and  dynamic  viscosity  in  the  Van  Driest  damping  factor 
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D  in  (13)  are  evaluated  at  the  wall,  rather  than  locally,  in  agreement  with 
most  studies  of  strong  shock-boundary  layer  interaction.^  ^  Second,  the 

g 

pressure  gradient  correction  to  N  has  been  omitted.  Previous  investigations 

of  flows  with  strong  shock-boundary  layer  interaction  have  been  performed 
13  10-12 

both  with  and  without  a  pressure  gradient  correction.  The  present 

code,  however,  employs  the  eddy  viscosity  relaxation  model  of  Shang  and 
Hankey^  (discussed  below) ,  which  has  been  developed  without  inclusion  of 
a  pressure  gradient  correction,  and  for  this  reason  the  pressure  gradient 
correction  has  been  omitted. 

The  transition  from  the  inner  equilibrium  eddy  viscosity  eeq.  to  the 
outer  eddy  viscosity  eeq  occurs  at  a  distance  at  which  eeq  =  eeq  . 

Thus,  the  full  form  of  the  equilibrium  eddy  viscosity  is 


eq 


eq, 


eq. 


n  <■  n 


n  >  n 


(15) 


The  transition  factor  F(s),  employed  to  allow  a  smooth  development  of 
the  eddy  viscosity  in  the  region  of  transition  from  laminar  to  turbulent 
flow,  is  given  by 


T(s)  =  1  -  exp(-0 .412s ^ ) 


s  =  (s  -  si)/A  , 


(16) 


Sj  <  s  <  Sj 


where  s^  and  s^  denote  initial  and  final  locations  of  the  transition 
zone  with 


A  =  slr=3/4  slr=i/4 
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As  applied  in  the  numerical  code,  the  transition  model  requires  only  the 

specification  of  and  s^.  These  values  may  be  obtained  by  the  method 
14 

of  Deem  and  Murphy  and  the  generally  accepted  criterion  that  the  Reynolds 

number  based  on  the  distance  from  the  leading  edge  approximately  doubles 

9 

across  the  transition  region. 

In  the  vicinity  of  the  interaction  of  a  shock  wave  with  a  turbulent 
boundary  layer,  a  modification  of  the  eddy  viscosity  is  required  in  order 
to  incorporate  the  experimental  observation  that  the  turbulent  shear 
stress  does  not  immediately  react  to  an  abrupt  change  in  the  flow.^  This 
effect  is  modeled  by  the  approach  of  Shang  and  Hankey,^’^  which  employs 
the  following  relaxation  model  for  the  eddy  viscosity: 

e(s,n)  =  eec](SQ,n)  +  [ce^(s,n)  -  (sq  ,  n)  ]  [  1  -  exp(-(s  -  Sq)/X)]  (17) 

where  s^  is  a  location  approximately  seven  boundary  layer  thicknesses 

upstream  of  the  intersection  of  the  shock  wave  and  the  wall,  c  is  the 

eq 

equilibrium  eddy  viscosity  given  by  (11)  to  (15)  and  X  is  typically  twenty 

times  the  boundary  layer  thickness  6  at  s  .  In  the  above,  s  is  the  distance 

along  the  surface  from  the  leading  edge,  and  n  is  the  distance  normal  to 

the  surt.ii'*.  The  above  expression  is  employed  for  s  >  s^,  and  provides 

for  a  relaxation  of  r  across  the  shock-boundary  layer  interaction  region, 

ultimately  approaching  the  equilibrium  expression  given  by  (11)  to  (15) 

for  s  -  s  much  greater  than  20  6  . 

o  o 

The  relaxation  model  is  important  at  shock-boundary  layer  interactions 
when  flow  separation  occurs.^’'*'  The  model  has  not  been  thoroughly 
investigated  for  boundary  layers  subject  to  repeated  shock  impingement 
(particularly  for  those  cases  where  the  distance  between  shock  impingement 
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is  less  than  the  relaxation  length  X),  and  has  therefore  been  incorporated 
in  the  numerical  code  for  use  at  one  shock-boundary  layer  interaction 


each  on  the  ramp  and  cowl. 

3.  Computational  Sublayer 

A  major  factor  governing  the  efficiency  of  the  numerical  algorithm 
employed  in  the  code  to  solve  the  Navier-Stokes  equations  is  the  requirement 
of  resolving  all  pertinent  scales  within  the  turbulent  boundary  layers  on 
the  ramp  and  cowl.  The  exceedingly  fine  mesh  spacing  required  to  resolve 
the  viscous  sublayer  portion  of  the  turbulent  boundary  layers  would 
ordinarily  severely  impair  the  efficiency  of  the  numerical  code.  In  order 
to  alleviate  this  difficulty,  a  separate  and  efficient  treatment  of  the 

k 

region  containing  the  viscous  sublayer  and  transition  portion  of  the 
boundary  layer  is  employed.  Within  this  region,  an  approximate  form  of  the 
full  equations  of  motion  is  utilized  in  agreement  with  previous  studies 
of  turbulent  boundary  layers. ^  The  governing  equations  are 


3u 1 

3y’ 


9JL_ 
3x ' 


3t 


3y’ 


x’y’ 


(18) 


7y~[™(cPT 


1  .2x 

2  u  )  - 


/  ^ 
Cp  Pr 


c 

Pr, 


,9^ 

3y' 


-  U  T 


=  0 


(19) 


x  y 


=  (u  + 


.  3u ' 

e)iy“ 


(20) 


where,  in  this  case,  x' 
normal  to  the  surface. 


and  y’  are  local  cartesian  coordinates 

respectively;  x  ,  ,  is  the  shear  stress 
x  y 


parallel  and 
parallel  to 


This  region  is  defined  by  0  <  n  <  50vw/u.,  where  v  =  u  / p  and 
u^  =  * Tw/Pw  •  In  this  context,  the  expression  "transition  region"  refers 
to  the  approximate  domain  lOv  /u.  <  n  <  50v  /u.  and  is  not  to  be  confused 
with  the  "transition  zone"  discussed  in  reference  to  (16). 
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the  wall;  and  u'  is  the  velocity  component  parallel  to  the  wall.  The 

equations  are  applied  within  a  narrow  region  0  £  y '  _<  adjacent  to  the 

walls  denoted  as  the  "computational  sublayer."  The  applicability  of  the 

above  expressions  (18)  and  (19)  imposes  an  upper  limit  on  the  value  of 

1 8 

y',  which  has  been  found  to  be  roughly  60v  /u.  .  The  achievement  of 
m  w  * 

maximum  code  efficiency  implies  that  the  value  of  y^  should  be  as  large 

as  possible  within  the  above  limit.  In  practice,  th ^  value  of  y^  is 

determined  by  the  mesh  obtained  from  the  coordinate  transformation 

program  and  the  nature  of  the  flowfield  (e.g.,  v  ,  t  ) ,  as  discussed 

w  w 

in  a  later  section.  Precise  control  of  y  at  every  location  is  not 

m  7 

afforded  by  the  numerical  codes.  However,  proper  operation  can  limit 
the  typical  range  of  y^  from  approximately  20\W u*  to  60v^/u*,  thereby 
achieving  an  efficient  and  accurate  solution  of  the  flow. 
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SECTION  II 


NUMERICAL  ALGORITHMS 
A.  Coordinate  Transformation 

The  function  of  the  coordinate  transformation  program  is  to  provide 
a  distributed  mesh  of  points  that  can  be  used  to  accurately  resolve  all 
pertinent  features  of  the  flow.  The  basic  concept  is  indicated  in  Figure 
3.  The  unit  square  in  the  transformed  plane  is  covered  by  a  uniform  mesh 
of  points  (c^,  n j )  where 


=  (i  -  1)  Ara  ,  i  =  1,  IL 

hj  =  (j  -  1)  An  ,  j  =  1,  JL 


(21) 


where  IL  and  JL  are  the  number  of  points  in  the  C-  and  n-directions, 

respectively.  By  definition,  Ac  =  1/ (IL  -  1)  and  An  =  1/(JL  -  1).  The 

values  of  Ac  and  An  are  constant,  although  not  necessarily  equal.  The 

image  of  the  points  (C^,  n^)  in  the  physical  plane,  obtained  from  equations 

(1)  and  (2),  is  a  general  curvilinear  mesh  of  points  with  corresponding 

cartesian  components  (x.  ,,  y.  .).  The  distribution  of  points  (x.  ,  y.  .) 

i,l  i,l  i,J  1,1 

is  highly  non-uniform  in  order  to  accurately  resolve  the  features  of  the 
flow.  The  governing  fluid  dynamic  equations  are  solved  in  the  transformed 
plane,  and  the  flowfield  within  the  inlet  is  obtained  directly  from  the 
knowledge  of  the  coordinate  transformation. 

In  order  to  facilitate  determination  of  the  coordinate  transformation, 
the  governing  equations  (1)  and  (2)  are  inverted.*  The  resultant  equations 
are 


12 


where 


32x 

26 

32x 

a2 

o  x  .  r 

3x 

o 

+  Y  — x  +  6 

ae2 

3?3q 

an2 

an 

.2 

_ 

26 

a2 

8  y 

+  Y  ^  +  6 

h. 

a?2 

3?3n 

3n2 

3n 

a 


2 


3x  3x  |  3y  3y 

3?  3n  3£  3n 


Y 


+ 


2 


(22) 


! 

I 

} 


(23) 


6  =  T(?,n)Y 

* 

where  T(5,rj)  is  given  in  (2) .  Equations  (22)  and  (23)  are  solved  in  the 
unit  square  in  the  transformed  plane.  The  coordinates  x(;,n)  and  y(c,n) 
are  subject  to  Dirichlet  boundary  conditions  on  the  boundaries  of  the 
unit  square,  i.e., 


x  =  f1(?) 

y  =  gj_U) 


on  n  =  0 


x  =  f2(c) 
y  =  g2(0 

x  =  f3(n) 
y  =  g3(n) 


on  n  =  1 


on  C  =  0 


•k 

The  symbols  u,  8,  Y  and  6  in  (23)  are  used  to  represent  different 
quantities  in  sections  other  than  the  present  (11. A). 
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x  =  f^(n) 


on  £  =  1 


y  =  8A(n) 


The  above  boundary  conditions  can  be  clearly  interpreted  as  specifying  the 
location  of  the  mesh  points  on  the  entire  boundary  of  the  computational 
domain  in  the  physical  plane  (e.g.,  the  ramp,  cowl,  upstream  and  downstream 
boundaries).  As  shall  be  discussed  later,  the  specification  of  the  boundary 
mesh  point  distribution  requires  careful  consideration  of  the  nature  of  the 
flowfield  to  be  computed. 

Equations  (22)  are  solved  numerically  using  an  Accelerated  Gauss-Seidel 

19  11 

Method.  Beginning  with  an  initial  estimate  (x.  . ,  y  .),  an  iterative 

l  , 1  1,3 

sequence  of  values  (x"  ,  y"  ),  n  =  2,3,  •••  is  obtained  from 

1,3  i  >  3 

n+1  n  ,  ~  n  . 

x,  ,  =  x  .  +  w,  . (x.  ,  -  x  ) 

1,3  1,3  1,3  1 » J  1,3 


n+1  n  .  ,  ~  n  . 

y.  .  =  y.  ,  +  w.  , (y.  .  -  y.  .) 

i,3  1,  j  i,j  i,3  i,3 


n  =  1,2,3, 


(24) 


where 


x  .  =  x?  .  +  a.  .(x?  .  .  +  xn+!  .  -  2x?  ,)/Ac2 

1,3  i,3  1,3  1+1, j  i-l,3  i,3 


_2r  (x11  -  xn  -  xn+1  +  x0^  )/4ArAn 

i,j  i+l,j+l  i-l,j+l  i+l,.j-l  i-lj-r7  fljAn 

(25) 

,  ,  n  ,  n+1  „  n  .  ,  2 

+  Y.  .(x.  . +  x.  .  ,  -  2x.  .)/An 
1,3  1,3+1  1,3-1  1,3 


+  <5.  Ax"  -  x1?"1"!  )/2An 

x,j  i , 3+1  i,3-l 


and  a  similar  expression  is  employed  for  . . 
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.1 


The  coefficients  a.  ,,  6.  . ,  y.  and  S.  .  are  obtained  by  straightforward 

i-  >  j  i ,  J 

second-order  accurate  finite-difference  expressions  for  the  terms  in  (23), 

e.g. , 


,,  n  n+1  ,2  ,  n  n+1  ,2  2 

a.  .  =  [(x.  -  x.  .  +  (v.  ...  -  v.  .  )  1/4  An 

i,J  i,j+l  i,J-l  i.j+l  i,l-l 

The  algorithm  (24)  is  used  to  update  the  coordinates  along  lines  of  constant 

p  (i.e.,  constant  j)  in  the  direction  of  increasing  r  (i.e.,  increasing  i). 

Thus,  in  forming  x.  .  and  y.  the  most  recent  values  of  x  and  v  on  the 
i-  *  J  i ,  J 

i  -  1  line  (e.g.,  x?"*".  ,,  yn+^  ,)  are  employed  to  facilitate  convergence 
i,j-l’  i,j-l 

and  reduce  computer  storage  requirements. 

The  acceleration  parameter  w .  .is  employed  to  speed  convergence  of 

i ,  J 

the  iteration  sequence.  Since  equations  (22)  are  non-linear,  an  optimum 
value  for  w  at  each  point  cannot  be  found.  A  linearized  stability  analysis 
of  the  highest-order  terms  yields  the  following  criterion  for  stability  of 
the  algorithm: 


.  „  2  2 

r\  ^  Ac.  an 

0  “l.j  <  2  "2 

An  «i,j  +  AA  A.j 

Based  on  experience  with  the  code,  the  following  criterion  is  employed: 

4  2  A  2 

w.  j  =0.95  - y— 

fin  ^  j  +  A(  j 

The  convergence  of  sequence  (24)  is  monitored  by  computing  the  following 


ratios  at  all  interior  points  at  selected  values  of  n: 
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(27) 


II 

•I— > 

•H 

i  n+1 
x.  . 

i  n+1 
x.  , 

1  iJ 

i.i 

i  n+1 
|yi,t 

-y"J/a  + 

,  n+1 

|yi,i 

These  quantities  essentially  represent  the  absolute  relative  change  in 
the  values  of  (x^  ^ ,  y^  during  the  previous  iteration,  with  an  additional 
additive  unitary  term  in  the  denominator  to  prevent  an  undefined  operation 


when  x.  .  or  y.  .  is  zero.  Iteration  sequence  (24)  is  considered  to  have 
i » J  i  9 1 

converged  when 


max  r  +  max  r  <  CONVER  (28) 

i,j  Xi,j  i >  j  yi,j 


where,  for  example,  max  rY  is  the  maximum  value  of  r„  over  all 

1»1  xi . j  xi,j 

interior  points  and  CONVER  is  a  small  quantity,  typically  10  or  smaller. 

Further  guidelines  for  the  choice  of  CONVER  are  given  in  Se'-'i'ion  III. C.  (Table  5). 

The  function  of  the  source  term  T(^,n)  in  (1)  is  to  provide  stretching 
of  the  mesh  in  the  p-direction  within  the  boundary  layers  on  the  ramp  and 
cowl  and  approximately  uniform  spacing  in  between  the  boundary  layers. 

The  effect  of  this  term  can  be  seen  by  considering  an  arbitrary  streamwise 
location  given  by  C  =  within  the  inlet.  For  typical  inlet  configurations, 
the  second  equation  in  (22)  may  be  approximated  to  lowest  order  by 


^y-  +  T«,p)^  =  0 

3n  3n 


(29) 


The  solution  to  this  equation  can  be  written  in  the  following  form,  noting 
that  Pj  =  (j  -  1 ) / ( JL  -  1): 
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yi  i  +  a1[exp(C1rij/Ti1)  -  1]  , 


(30) 


y 


i,  j 


<  1  +  a^fexpCC^)  -  1]  +  a9(n  -  n  j )  ,  ''l  —  n  -  n2 

yi  jl  -  a3texp(a4(l  -  ru))  -  1],  n9  '  n.  fl 


(31) 

(32) 


where 


_  _L  ,  vJ  _L  ,  -L 

ao  n1[exp(C1)  -  1 ] /C^  +  exp(C1)[n2  -  +  (1  -  n9)(l  -  exp(-C2))/C9] 

al  =  aoVCl 

a2  =  aQ  expCC^)  (33) 

a3  =  a0(exP(ci  _  C2))(l  -  P 2 ) / C 9 
a4  =  C27 ^  _  n2^ 

Noting  that  the  expressions  a^,  ,  a2>  a^,  and  terms  C^,  0 9 ,  and  n2 

are  not  functions  of  n,  it  is  evident  that  the  above  approximate  solution 

for  y.  .  is  characterized  by  an  exponentially  stretched  mesh  between 
1 ,  J 

n  =  0  and  n  =  n^,  a  uniformly  distributed  mesh  between  n  -  and  n  =  n9, 
and  an  exponentially  stretched  mesh  between  n  =  1  and  n  =  n9  <  1. 

The  application  of  the  source  term  T(5,n)  requires  specification  of 
the  quantities  ,  C^,  and  n2  at  each  value  of  r  ,  i  =  1,2,  IL. 

In  practice,  this  is  accomplished  in  the  following  fashion.  First,  values 
for  and  n2  are  chosen  that  reflect  the  desired  number  of  points  to  be 
stretched  in  the  vicinity  of  the  upper  and  lower  surfaces.  By  definition, 

=  (j  -  1)/(JL  -  1),  where  j  is  the  number  of  exponentially  stretched 
points  near  the  lower  boundary  (including  the  point  on  the  boundary),  and 
n2  =  (.jj  -  1)/(JL  -  1),  where  the  number  of  exponentially  stretched  points 


near  the  upper  boundary  is  JL  -  j^  +  1.  In  conjunction  with  a  separate 

program  (called  BNDRY) ,  which  is  discussed  in  a  subsequent  section,  the 

values  of  1  and  1,  are  chosen  in  order  to  provide  a  sufficient  number  of 
Jo  1 

points  within  the  boundary  layers  on  the  ramp  and  cowl.  Although  the 
coordinate  transformation  program  permits  j  and  j^  to  rary  with 
satisfactory  results  have  been  obtained  with  judiciously  chosen  constant 
values  for  j  and  j,. 

Secondly,  the  values  of  C^(C)  and  02(c)  are  determined  by  the  requirement 
that  the  height  above  the  wall  of  the  first  row  of  points  adjacent  to  the 
upper  and  lower  surfaces  (i.e.,  q  =  An  and  q  =  1  -  An ,  which  correspond 
to  j  =  2  and  j  =  JL  -  1 ,  respectively)  remains  within  the  limits  indicated 
by  the  computational  sublayer  model  (see  Section  III. A)  at  all  streamwise 
locations  C .  This  is  accomplished  automatically  by  the  coordinate  trans¬ 
formation  code,  using  information  specified  by  the  user  that  relates  the 
approximate  desired  height  of  the  rows  q  =  An  and  n  =  1  -  Aq  from  the 
lower  and  upper  surfaces,  respectively,  as  a  function  of  streamwise  position 
C  (see  Section  III)  .  From  (30)  the  approximate  height  of  the  first  row  of 
pc ints  above  the  lower  surface  is  the  following: 

yi  2  "  yi,l  =  a1IexP(c1  An/h!)  “  1]  (34) 

while  the  height  of  the  row  of  points  adjacent  to  the  upper  boundary  is 

yi , JL  ~  yi,J2  =  a3f  exP(a4  Aq)  -  1)  (35) 

The  above  expressions  are  non-linear  equations  for  C-j  and  C^,  and  are  solved 
19 

by  Newton's  method. 
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In  solving  the  governing  fluid  dynamic  equations  (3)  to  (6),  the 

coordinate  transformation  derivatives  5  ,  5  ,  n  .  n  are  required.  They 

x’  y  x  y 

may  be  expressed  as 

3C  =  _i_  3Z 
3x  J  ’  Drj 

3  5 _ 1  3x 

3y  J’  3q 

(36) 

3n  =  _  _L  ii 

3x  J’  dr, 

3n  _  1  3x 
3y  J'  dr, 

where  J'  =  (3x/3£) (3y/3n)  -  (3x/3q) (3y/30  -  At  all  interior  points,  the 
derivatives  on  the  right-hand  side  of  (36)  are  obtained  using  second-order 
central  differences,  while  on  the  boundaries  of  the  domain  second-order 
one-sided  differences  are  employed  as  required.  In  addition,  the  coordinate 
transformation  code  permits  the  recognition  of  a  finite  number  of 
discontinuities  in  the  slope  of  the  upper  and  lower  surfaces  between 
£  =  0  and  £=1,  introducing  additional  second-order  one-sided  differences 
to  account  for  them. 

B.  Navier-Stokes  Code 
1 .  MacCormack's  Method 

The  inlet  flowfield  is  obtained  by  numerical  integration  of  the 

Navier-Stokes  equations  (3)  in  time  from  an  assumed  initial  condition  until 

a  steady-state  solution  is  obtained.  The  numerical  algorithm  employed  is 

20  21 

MacCormack's  method,  ’  which  is  an  alternating-direction  explicit 
technique  of  Lax-Wendroff  type.  It  has  been  applied  to  the  solution  of  a 
wide  variety  of  problems  in  high  speed  flow  involving  strong  viscous-inviscid 
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r 


interactions.  A  representative  sample  of  applications  is  indicated  in 

References  10,  16  and  20  through  30. 

The  numerical  algorithm  is  applied  on  the  finite-difference  mesh  in 

the  U ,n)  transformed  plane  indicated  in  Figure  3.  Denoting  a  sequence 

of  time  levels  tn  =  n  At,  n  =  1,2,3,  ••*,  where  At  is  the  time  step,  the 

vector  of  dependent  variables  U?  .  at  position  (t.,n.)  and  time  tn  is 

1,3  i  1 

updated  to  time  tn+^  by  the  following  expression: 

Un+1  =  L(At)Un  . 

i,l  i,3 

where  L( At)  is  a  symmetric  sequence  (discussed  below)  of  time-split 
one-dimensional  difference  operators  L  (At  )  and  L  (At  ) .  The  operator 

C  r,  n  n 

is  a  second-order  accurate  finite-difference  algorithm  for  the 
one-dimensional  equation 


1W+1£  =  o 

3t  H 


*  ** 

Using  the  dummy  time  indices  and  with 


**  * 

U.  =  L  (At  )U.  . 
i,j  C  ?  i,J 

the  L operator  is  given  by  the  following  two-step  predictor-corrector 
me  thod : 


**  *  a  r  •k  * 

Predictor  Step:  U.  .  =  U.  ,  -  — — ^(F.  ,  -  F.  .  .) 

i,J  i,j  i,j  l-l, j 


**  I  *  **  “  r  **  ** 

Corrector  Step:  U.  .  =  —  U.  .  +  U.  .  - ^(F...  .  -  F.  .) 

i,l  2  1,1  1,1  Ar.  i+l,i  i,3 


*  * 

where  r  implies  that  the  flux  terms  are  evaluated  using  U.  .  and  so 
i,j  i,3 

forth.  Further  information  regarding  the  differencing  of  the  stress  terms 


in  F  is  given  in  References  21  and  25.  The  operator  L  (At  )  is  a  second- 

n  n 

order  finite-difference  algorithm  for  the  one-dimensional  equation 


St  3n 


Using  a  similar  dummy  index  notation  as  in  (38),  the  L  operator  is 
given  by  the  following  two-step  predictor-corrector  method: 

Predictor  Step:  (J .  .  =  U  -  - — — (G.  .  -  G.  ) 

1,3  i,j  An  1,3  1,3-1 


**  1 

Corrector  Step:  U.  .  =  -r- 
i>3  2 


** 


** 


At 


** 


kk 


U.  .  +  U.  .  -— '-(G.  -  G.  .) 

1,3  1,3  An  1,3+1  1,3 


where,  for  example,  G.  implies  that  the  flux  term  is  evaluated  using 

_  1,3 

** 

U .  . ,  etc . 

i,j 

The  operator  l(At)  in  (37)  is  constructed  from  a  symmetric  sequence 

of  the  operators  L  (At)  and  L  (At  ).  In  the  numerical  code,  two 

C  C  n  n 

particular  forms  are  employed,  specifically 


L(At)  =  L  (At/ 2) L  (At) L  (At/2) 

n  An 


(401 


and 


L(At)  =  L  (At/2) L  (At) L  (At/2) 
t  nr. 


(41) 


where  the  operators  on  the  right-hand  side  are  applied  to  (/ .  ,  sequentially 

from  right  to  left.  Both  sequences  (40)  and  (41)  provide  a  spatially  and 
temporally  second-order  accurate  integration  of  the  ful  !  Navi er-S tokos 
equations  (3).  The  time  step  size  indicated  hv  the  argument  of  each 
operator  L ^  and  must  not  exceed  the  maximum  allowed  for  that  operator. 

An  approximate  linearized  stability  analysis  yields  the  following,: 
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(42) 


At 


C 


min  At. 


At 

n 


min  At. 


(4  3) 


where 


As 

_ 5. 


|u  1  +  c  + 

261  ^2_" 

/p 

1  C1 

As 

As 

L  C 

nJ 

_ 

(44) 


As 


At. 


— 

r  26, 

e„ ' 

- 

1  u  +  c  + 

1  n 

1 

As 

L  n 

/p 

(45) 


i.j 


where  min  denotes  the  minimum  value  of  the  right-hand  side  evaluated  at 

i,  j 

all  points  (t.,n,)  within  a  given  region  of  the  flow,  and 
■*-  J 

/~ ?  T 

As?  =  A  +  Cy 


As 


=  An//2  ■  2 


n  +  n 
x  y 


u  =  (uC  +  VC  )/  //  +  C2 
C  x  y  x  y 


(46) 


u  =  (un  +  vn  )/ /i2  +  n 2 
n  x  y  x  y 


c  =  /yRT  , 


R  =  c  -  c  =  gas  constant 
p  v 


8^  -  max 

1 2(p  +  e)  +  XT|  , 

+  M\ 

_ 

l  t>A 

9.  =  /| (p  +  e)A 
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In  practice,  the  quantities  (44)  ,  (45)  are  multiplied  by  a  factor  (denoted 
as  CFL)  that  is  slightly  less  than  one. 

2 .  Mesh  Splitting 

Due  to  the  non-uniformity  of  the  mesh  spacing  in  the  physical  plane, 
the  values  of  Ati  and  At?  .  vary  substantially  within  the  flow. 

Xi,j 

Specifically,  the  small  mesh  spacing  in  the  n-direction  within  the  boundary 
layers  (which  implies  large  values  of  and  consequently  small  values  of 
As  )  results  in  particularly  small  values  of  At-?  .  On  the  other  hand, 

n  i,j 

the  larger  values  of  As^  occurring  outside  the  boundary  layer  imply  larger 
values  of  At2.  ..  If  the  min  in  (43)  were  taken  over  all  points  in  the 

1,J  i,j 

flow,  the  value  of  At  would  be  controlled  by  the  values  of  Ato  in  the 

n  2i,.1 

boundary  layer,  and  hence  for  a  large  number  of  points  outside  the  boundary 
layers  the  operator  in  (40)  or  (41)  would  be  applied  using  a  value  of 
At^  that  is  smaller  than  the  maximum  allowed  locally  by  (45).  This  implies 
that  the  algorithm  is  not  being  optimally  applied. 


A  substantial  improvement  is  realized  by  incorporation  of  the  split 
21 

mesh  technique.  As  utilized  in  the  numerical  code,  the  computational 
domain  is  divided  in  the  q-direction  into  five  regions  (see  Figure  4) 
given  by 


Region 
1 
2 

3 

4 

5 

where  1  <  JI1  <  JI2  <  JI3  <  JI4  <  JL . 
sequence  is 


Extent 
1  <  j  <  JI1 
JI1  +  1  <  j  <  .112 
JI2  +  1  <  j  <  .113 
JI3  +  1  <  j  <  JT4 
JI4  +  1  ^  j  5  .1L 

Within  each  region  the  operator 
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,  £  =  1,2,4, 
(47) 


Regions  1,2,4  and  5: 


l(At)  =  [L  (At/ 2m  )L  (At/m.) L  (At/2m  ) ] 

n  £  4  x.  n  x. 


Region  3:  L ( A t)  =  (At/2) (At) L  (At/2)  (48) 

where  m  ,  £  =  1,  2,  4,  5  are  Integers,  and  the  exponent  m^  implies  that 
the  operator  sequence  within  the  brackets  is  applied  m^  times.  The 
application  of  the  stability  requirements  (42)  to  (45)  within  each  region 
yields  the  following: 

Regions  1,  2,  4,  5:  At/m^  £  min^At^  ’  ^  (49) 

i,j  i  >  j 

Region  3:  At  £  min^^At^  ,  At^  )  (50) 


where  min„(At,  ,  2At0  )  implies  the  minimum  value  of  At,  and 

X,  -t-4  A  ^  •  4  i,'  4 


i.j 


2Ato  within  the  £  region.  A  further  requirement  is  that  the  ratios 
rn^/n^  and  m^/m^  be  integers,  in  order  that  adjacent  regions  may  be 
updated  in  time  in  a  near-simultaneous  fashion. 

An  automatic  optimization  algorithm  is  utilized  in  order  to  determine 
the  most  efficient  splitting  of  the  five  regions.  At  each  time  step, 
the  following  ratio  is  computed  for  every  set  of  values  of  JI1,  JI2, 

JI3  and  JI4: 


At 


5 

l 

t  =1 


P„m„ 
£  £ 


(51) 


where  is  the  number  of  rows  in  region  £  (e.g.,  =  JI2  -  JI1)  and 

is  defined  in  (47)  with  m^  =  1 .  This  quantity  is  essentially 


proportional  to  the  ratio  of  the  time  step  At  to  the  computer  time  required 
to  update  the  flowfield  by  the  amount  At,  and  hence  is  a  direct  measure  of 
the  relative  efficiency  of  a  particular  mesh  splitting.  The  set  of  values 
of  JI1,  JI2,  JI3  and  JI4  that  yields  the  largest  value  of  the  ratio  (51) 
while  having  integer  values  for  m^/rr^  and  m,_/m^  is  then  employed  in 
updating  the  solution  by  the  amount  At.  The  process  is  continually 
repeated,  thereby  achieving  maximum  code  efficiency  at  all  times.  To 
the  author's  knowledge,  this  represents  the  first  incorporation  of  a 
completely  automatic  mesh-splitting  algorithm  with  MacCormack's  method. 

Special  consideration  is  required  at  the  intermesh  boundaries  between 
the  five  regions  in  order  to  conserve  mass,  momentum  and  energy.  Details 
of  the  intermesh  conservation  relations  are  given  in  Reference  21. 

3 .  Mesh  Overlapping 

Due  to  limitations  in  computer  core  memory,  it  is  often  not  feasible 
to  compute  the  entire  flowfield  of  a  high  speed  inlet  using  a  single  mesh. 
In  such  cases,  the  technique  of  mesh  overlapping  is  useful.  The  basic 
concept  is  illustrated  in  Figure  5.  The  entire  inlet  flowfield  is  divided 
into  two  or  more  overlapping  regions  denoted  by  A,  B,  etc.  A  steady-state 
solution  is  obtained  first  in  region  A  (defined  by  the  boundary  abed)  .  The 
values  of  the  flowfield  variables  at  station  a'd'  are  then  used  as  the 
upstream  profile  for  region  B  (defined  by  the  boundary  a'b'c'd').  The 
equations  of  motion  are  then  integrated  in  time  in  region  B  until  a 
steady-state  solution  is  achieved.  Further  downstream  regions  can  be 
computed  in  the  same  fashion. 

The  application  of  the  mesh  overlapping  technique  requires  that  two 
conditions  be  satisfied.  First,  within  each  overlapping  zone  (e.g.. 


1* 


1 


a’b  c  d'),  the  flowfield  outside  the  boundary  layers  must  be  supersonic 
in  order  that  there  be  no  downstream  influence  on  the  flow  in  this  area. 
Secondly,  within  the  vicinity  of  the  restart  station  (curve  a'd’),  the 
boundary  layers  on  the  ramp  and  cowl  must  be  developing  smoothly.  In 
particular,  the  restart  station  must  not  coincide  with  any  strong 
viscous-inviscid  interaction  (e.g.,  shock-boundary  layer  interaction)  or 
abrupt  change  in  conditions  in  the  boundary  layer  (e.g.,  an  abrupt  change 
in  the  bleed  mass  flux) .  This  condition  implies  that  the  behavior  of  the 
flow  in  the  boundary  layer  in  the  vicinity  of  the  restart  station  is 
essentially  governed  by  the  boundary  layer  equations  with  no  strong 
viscous-inviscid  coupling,  and  hence  mathematically  there  is  no  downstream 
influence  on  the  flow  in  the  boundary  layers  near  the  restart  station. 

The  technique  of  overlapping  mesh  regions  has  been  widely  used  in 
high  speed  flow  calculations,  including  inlet  f lows  ^ ^ ^  ^as 

been  tested  and  proven  in  the  present  numerical  code. 

4 .  Boundary  Conditions 

The  boundary  conditions  for  the  Navier-S tokes  equations  can  be 
categorized  into  four  major  types  (see  Figure  6a) — upstream,  downstream, 
wall  and  no  reflection.  On  the  upstream  boundary,  the  flow  variables  are 
held  fixed  at  the  appropriate  freestream  values.  In  the  case  of  mesh 
overlap,  the  variables  at  the  restart  station  (a'd'  in  Figure  5)  are 
held  fixed  at  their  steady-state  values.  At  the  downstream  boundary,  the 
conventional  zero-gradient  boundary  condition  is  applied,  i.e.. 


On  the  ramp  and  inlet  surfaces,  the  following  boundary  conditions  are  used 


w 


v  *  s  =  0  (Zero  tangential  velocity 

— ►  -V 

v  •  n  =  (Boundary  layer  bleed) 


3T 

3n 


=  0  (Adiabatic  wall) 


3n 


=  0  (Approximate  derived  boundary  condition 
for  pressure) 


(53) 


where  s  and  n  are  tangential  and  normal  unit  vectors,  respectively,  at 

the  wall  that  can  be  expressed  in  terms  of  n  and  n  .  Also,  v  is  the 

x  y  w 

normal  velocity  at  the  wall,  obtained  from  the  specified  mass  flux 

along  the  boundaries.  The  approximate  derived  boundary  condition  for 

the  pressure  has  been  successfully  employed  in  a  variety  of  flows 
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exhibiting  strong  viscous-inviscid  interaction.  ’  ’  The  fourth 

category  of  boundary  conditions  refers  to  curves  AB  and  DE  of  Figure 

6a.  As  the  name  suggests,  these  contours  are  assumed  to  be  no-reflection 

30 

boundaries.  The  boundarv  conditions  are 


3u 

3£ 


=  0 


3v 

3S 


=  0 


3T 

3? 


=  0 


i£ 


0 


(54a) 


where  the  derivative  3/  3£;  is  taken  along  the  outwards  running  characteristic 
at  the  boundary,  which  is  oriented  at  the  Mach  angle  n  =  sin  ^(1/M)  with  respect 
to  the  velocity  vector  v  as  illustrated  in  Figure  6b.  The  velocity  v  may, 
in  general,  be  pointing  in  or  out  of  the  computational  domain  at  either 
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boundary.  The  utilization  of  the  no-reflection  boundary  condition  requires 
that  the  following  conditions  be  satisfied: 

a)  The  local  Mach  number  M  =  jv|/c,  where  c  is  the  local  speed  of 
sound,  must  be  greater  than  one. 

b)  The  local  normal  Mach  number  |v  •  n|/c  must  be  less  than  one. 

The  first  requirement  implies  that  the  characteristics  exist  (i.e.,  the 
flow  is  locally  supersonic) ,  and  the  second  requirement  insures  that  there 
is  a  single  outwards  running  characteristic.  The  no-reflection  boundary 
condition  permits  the  angle  of  attack  of  the  inlet  to  be  varied  by 
constructing  the  contours  AB  and  DE  to  be  parallel  to  the  incoming  flow, 
with  the  cartesian  x-axis  aligned  with  AB  and  DE .  In  those  instances  where 
the  shock  wave  emanating  from  the  ramp  leading  edge  E  intersects  the  lower 
surface  upstream  of  the  cowl  leading  edge  B,  the  contour  between  the 
intersection  point  and  B  may  have  to  be  roughly  aligned  with  the  local 
flow  direction  in  order  to  satisfy  condition  b)  above.  If  either  condition 
a)  or  b)  is  not  satisfied  locally,  the  numerical  code  assigns  freestreara 
values  to  the  flow  variables  according  to 


u  =  U 

GO 


T  =  T 

CO 


and  prints  a  warning  message.  Since  the  boundary  conditions  (54b)  may 
not  be  physically  realistic  under  such  circumstances,  the  user  must  take 
appropriate  steps  to  remedy  the  situation  if  the  warning  messages  persist 
throughout  the  calculation.  If  condition  a)  is  being  violated  (e.g.,  a 


terminal  shock  is  located  upstream  of  the  cowl  lip,  as  occurs  in  fully 
external  compression  inlets),  it  should  be  recognized  that  the  boundary 
conditions  (54)  are  inappropriate  and  the  flowfield  cannot  be  computed 
unless  proper  modifications  are  made  to  the  boundary  condition  subroutine. 

If  condition  b)  is  not  satisfied  (e.g.,  the  ramp  shock  intersects  the 
curve  AB  upstream  of  the  cowl  leading  edge,  and  the  normal  Mach  number  on 
the  segment  between  the  intersection  point  and  cowl  leading  edge  is  greater 
than  one) ,  the  user  can  realign  the  appropriate  segment  of  the  boundary  as 
discussed  previously. 

The  boundary  conditions  are  implemented  in  the  numerical  code  using 
second-order  accurate  finite-difference  approximations  in  general .  The 
zero  normal  derivative  boundary  condition  on  the  pressure  in  (53)  is 
incorporated  using  first-order  accurate  differencing. 

5 .  Numerical  Damping 

The  algorithm  of  MacCormack  is  a  "shock-capturing"  type  in  which 
shock  waves  are  effectively  broadened  or  "smeared"  over  several  mesh  points, 
rather  than  appearing  as  sharp  discontinuities  in  the  flow.  In  the  vicinity 
of  strong  shock  waves,  the  numerical  truncation  error  associated  with  the 
shock  broadening  can  lead  to  numerical  instability.  Although  this  difficulty 
can,  in  principle,  be  avoided  by  refining  the  mesh  sufficiently  in  the 


vicinity  of  the  shock  to  resolve  its  physical  structure  (i.e.,  by  employing 

mesh  spacing  normal  to  the  shock  of  the  order  of  the  mean  free  path  of  the 

gas),  the  cost  would  be  prohibitive.  Alternately,  a  fourth-order  "pressure 

21 

damping"  term  is  employed,  which  is  of  significance  only  in  the  vicinity 
of  pressure  oscillations  where  the  solution  has  been  adversely  affected  by 


truncation  error  anyway. 


In  the  L 

4 


operator  in  (40)  and  (41),  the  following 
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term  is  added  to  the  flux  F.  . 


•k-k 

where  the  expression  is  evaluated  using  level  quantities.  Corresponding 

expressions  are  used  for  the  1^  operator.  The  numerical  code  permits  the 

user  to  omit  the  pressure  damping  terms  for  the  purpose  of  increased 

efficiency  when  computing  flows  with  no  shock  waves. 

An  additional  form  of  damping  is  incorporated  into  the  numerical  code 

in  order  to  control  any  destabilizing  oscillations  that  may  result  from 

32 

initial  flow  transients.  The  technique  consists  of  multiplying  the 
viscosity  =  — (2/3) (u  +  e)  by  a  large  negative  constant  8  during  the 
early  stages  of  the  calculation,  resetting  A^,  to  -(2/3) (u  +  c)  (i.e., 

8  =  +1)  and  then  converging  the  flow.  In  all  calculations  with  the 
numerical  code,  it  was  not  found  necessary  to  utilize  this  damping. 

Finally,  the  convective  damping  technique  of  MacCormack ^ is 
incorporated  into  the  operator  in  order  to  prevent  occurrence  of  a 
non-linear  instability  in  regions  of  separated  flow. 
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6 .  Convergence  of  Flowfield  to  Steady  State 

Since  the  steady-state  flowfield  is  obtained  by  integration  in  time 
from  an  assumed  initial  condition,  a  criterion  for  convergence  must  be 
stipulated.  Appealing  to  the  physics  of  the  flow,  it  is  evident  that 
the  physical  time  (referring  to  the  time  t  in  equation  (3))  necessary 
for  achievement  of  a  steady-state  solution  is  a  multiple  of  the  character¬ 
istic  time  t  required  for  a  fluid  particle  to  ♦ravel,  from  the  upstream 
to  the  downstream  boundaries  of  the  computational  domain.  Denoting  an 
approximate  value  of  the  velocity  in  the  inviscid  region  by  and  the 
length  of  a  given  region  by  L, 

t  =  L/U  (57) 

c  e 

Experience  has  shown  that  a  total  physical  time  of  approximately  3t  is 

sufficient  to  guarantee  achievement  of  a  steady-state  solution  provided 
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there  is  no  flow  separation.  For  cases  with  flow  separation,  a  somewhat 

larger  physical  time  of  5tc  to  6t  is  required.  The  above  criteria  are 

general  guidelines.  In  application,  a  prudent  strategy  is  to  check  the 

change  in  the  flow  variables  (e.g.,  surface  pressure,  skin  friction) 

between  two  physical  time  levels  that  are  separated  by  approximately  It' 

(e.g.,  t  =  2t  and  t  *  3t  ) .  Provided  that  the  changes  are  tvpicallv 
c  c  * 

one  to  two  percent  or  less,  the  flowfield  may  be  considered  converged. 

7 .  Computational  Sublayer 

The  geometry  associated  with  the  computational  sublayer  calculation 
is  illustrated  in  Figure  7  for  the  n  =  0  and  q  =  1  surfaces.  The  "ordinary 
mesh"  obtained  from  the  coordinate  generation  prot  am  is  indicated  bv  the 
open  symbols.  The  algorithm  of  MacCormack  is  employed  on  the  ordinary 
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mesh.  The  region  between  the  wall  and  the  first  row  of  adjacent  ordinary 
mesh  points  is  denoted  as  the  "computational  sublayer."  Within  this 
region,  a  locally  orthogonal  coordinate  system  (x',y')  is  employed.  Tire 
approximate  governing  equations  (18)  to  (20)  are  solved  on  a  finite-difference 
mesh,  indicated  by  the  closed  symbols. 

As  indicated  previously,  the  achievement  of  maximum  code  efficiency 

implies  that  the  distance  y^  of  the  first  row  of  ordinary  points  should  be 

as  large  as  possible  within  the  upper  limit  of  roughly  60 v  /u  .  Since 

w  * 

this  height  is  substantially  greater  than  the  viscous  sublayer,  accurate 
values  for  the  stress  components  and  heat  transfer  vector  cannot  be 
obtained  at  the  matching  points  by  ordinary  finite  differences  involving 
the  velocity  components  at  the  wall  and  the  matching  points.  The  purpose 
of  the  computational  sublayer  technique  is  to  provide  values  of  the  stress 
components  and  heat  transfer  vector  at  the  wall  and  at  the  matching  points 
for  use  in  application  of  MacCormack's  method  at  the  matching  points. 

The  cartesian  stress  components  x  ,  x  ,  x  are  obtained  from  the 
sublayer  solution  by  a  simple  coordinate  transformation,  under  the 
reasonable  assumption  that  the  normal  stresses  in  the  x'-y’  coordinate 


system  (i.e.,  Tx'x»’Tyiy<)  are  negligible  compared  to  Tx?v>-  Similarly, 
the  heat  transfer  components  Q^,  are  obtained  by  coordinate  transformation 
from  the  sublayer  solution  under  the  reasonable  assumption  that  <<Q^, 
in  the  sublayer. 


The  computational  sublayer  region  is  solved  after  each  updating  of 


the  adjacent  region  of  the  ordinary  mesh.  Denoting  Txiyi  as  T»  with 
subscripts  w  and  m  implying  the  wall  and  matching  point,  respectively,  the 


sequence  of  solution  is  as  follows.  First,  the  shear  stress  components 
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are  found  from 
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where  u^  is  the  component  of  the  velocity  at  the  matching  point  parallel 
to  the  surface  and 


f  (y  ’) 


exp 


m  dy  ’ 
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y  ’ 


(59) 


Expression  (58)  is  obtained  by  integrating  (18)  twice.  The  velocity 
in  the  sublayer  is  then  obtained  by  integrating 


3u  ’  _  m 
3y 1  (p  +  e) 


'(»'  +  e) 


(60) 


The  temperature  field  is  then  determined  by  (19),  and  the  heat  transfer 
vectors  computed.  Further  details  are  given  in  Reference  18. 

C.  Overall  Description  of  Programs 

There  are,  in  general,  a  total  of  four  separate  computer  programs 
that  are  employed  to  compute  the  flowfield  In  a  high  speed  inlet.  The 
purposes  of  the  programs  are  as  follows: 

a)  Coordinate  Generation  Program  (COORD) 

This  program  is  employed  to  generate  the  sequence  of  overlapping 
mesh  systems  that  are  utilized  bv  the  Navi er-Stokes  code. 


nm 


b)  Upstream/Downstream  Boundary  Mesh  Distribution  Program  (BNDRY) 

The  purpose  of  this  code  is  to  determine  the  appropriate 
distribution  of  mesh  points  in  general  for  the  upstream  and 
downstream  boundaries  of  each  mesh  system.  It  is  utilized  in 
conjunction  with  the  coordinate  generation  program. 

c)  Navier-Stokes  Program  (INLET) 

This  program  solves  for  the  steady-state  flowfield  within 
each  mesh  system. 

d)  Upstream  Boundary  Interpolation  Program  (UPSTRM) 

The  purpose  of  this  short  program  is  to  permit  interpolation 

of  flow  data  at  the  restart  station  of  two  overlapped  meshes 

where  the  height  of  the  two  mesh  systems  is  not  identical.  This 

procedure  is  employed  in  the  calculation  of  the  simulated  high 
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speed  inlets  developed  by  McDonnell  Aircraft  Company  (MCAIR) 
and  is  illustrated  in  Figure  8. 

For  the  MCAIR  inlet  configuration,  the  calculation  proceeds  in  two 
basic  stages.  In  the  first  stage,  the  boundary  layer  development  on  the 
ramp  (i.e.,  the  lower  surface  in  Figure  8)  is  computed.  This  procedure 
is  accomplished  in  the  following  steps. 

a)  Using  program  BNDRY,  the  desired  mesh  spacing  on  the  upstream  and 
downstream  boundaries  of  mesh  A  is  determined. 

b)  Using  the  information  from  a)  above,  the  coordinate  system  of 
mesh  A  is  computed  using  program  COORD. 

c)  Using  the  coordinate  transformation  computed  in  b) ,  the  steady-state 
flow  in  region  A  is  computed  using  program  INI, FT. 
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In  the  event  that  the  entire  ramp  boundary  layer  development  upstream  of 
the  inlet  entrance  cannot  be  computed  with  a  single  mesh  system,  a  second 
overlapping  mesh  system  is  generated  using  COORD  and  the  flowfield  is 
determined  using  INLET.  This  procedure  is  repeated  until  the  ramp  boundary 
layer  has  been  computed  up  to  the  vicinity  of  the  inlet  entrance. 

In  the  second  stage,  the  flowfield  within  the  inlet  is  computed. 

This  procedure  requires  the  following  steps: 

a)  Using  program  BNDRY ,  the  desired  mesh  spacing  on  the  upstream  and 
downstream  boundaries  of  mesh  B  is  determined.  See  also  step  c) 
immediately  below. 

b)  Using  this  information,  the  coordinate  system  of  mesh  B  is  determined 
using  program  COORD. 

c)  Because  the  vertical  distributions  of  points  in  meshes  A  and  B 
at  the  restart  station  are  not  identical,  the  program  1TSTR11  is 
used  to  interpolate  the  flow  variables  onto  the  upstream  boundary 
of  mesh  B.  There  are  certain  restrictions  on  the  use  of  program 
UPSTRM  (see  Section  IV.u.l)  that  must  be  satisfied  in  determining  the 
mesh  distribution  on  the  upstream  boundary  of  mesh  B. 

d)  The  steady-state  flowfield  is  obtained  using  program  INLET. 

For  all  subsequent  regions,  the  following  steps  are  taken: 

a)  The  program  BNDRY  is  used  to  determine  the  desired  mesh  spacing 
on  the  downstream  boundary  on  the  next  mesh  region. 

b)  The  coordinate  transformation  is  obtained  using  COORD. 

c)  The  steady-state  flowfield  is  obtained  using  INLET. 

In  those  cases  where  a  variety  of  inlet  configurations  are  to  be  tested  at  the 
same  entrance  conditions,  the  first  stage  need  only  he  accomplished  once. 
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SECTION  III 


GENERATION  OF  COORDINATE  TRANSFORMATION 
A.  Preliminary  Analysis 

The  basic  requirement  for  a  f inite-dif ferenci  mesh  is  to  provide 
accurate  resolution  of  all  pertinent  features  in  the  flow.  The  mesh 
point  distribution,  therefore,  is  dependent  on  the  particular  flowfield 
to  be  computed,  and  requires  preliminary  analysis  in  two  major  areas, 
as  discussed  below. 

1 .  Streamwise  Mesh  Spacing 

In  general,  the  controlling  influence  on  streamwise  mesh  spacing  for 

high  speed  inlet  calculations  is  the  necessity  of  accurately  resolving  the 

boundary  layer  development.  Unfortunately,  there  does  not  exi? t  precise 

criteria  for  the  determination  of  the  appropriate  streamwise  mesh  spacing, 

except  an  exhaustive  (and  usually  infeasible)  truncation  error  study  in 

each  case.  An  oft-quoted  guideline  in  computational  fluid  mechanics  is 

that  a  streamwise  mesh  spacing  with  cell  Reynolds  number  less  than  two 

is  sufficient  to  resolve  all  characteristic  streamwise  features  within 
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the  boundary  layer.  However,  this  criterion  is  based  on  a  truncation 
error  analysis  for  the  particular  case  in  which  first-order  upwind 
differencing  is  employed  for  the  convective  terms  of  the  governing  equations, 
and  does  not  specifically  apply  to  MacCormack 's  method.  Furthermore, 
experience  in  strong  viscous-inviscid  interacting  boundary  laver  flows 
has  shown  this  criterion  to  be  unnecessarily  severe/''"’^ 

A  generally  successful  guideline  for  determining  st reamwise  mesh 
spacing  is  based  on  the  ratio  Ax/4,  where  Ax  is  the  streamwise.  mesh  spacing 
(typically  evaluated  on  the  upper  and  lower  boundaries)  and  4  is  a  reference 
boundary  layer  thickness.  In  those  regions  in  which  the  boundary  is 


approximately  in  equilibrium  (e.g.,  mesh  A  of  Figure  8),  a  mesh  spacing 

of  Ax/6q  :  1  to  2  is  satisfactory,  where  6q  is  the  boundary  layer  thickness 

at,  say,  the  downstream  boundary  of  the  mesh  region.  For  regions  of 

strong  viscous-inviscid  interaction,  such  as  occur  within  the  inlet,  the 

value  of  Ax/6q  (where  6o  is  the  boundary  layer  thickness  immediately  upstream 

of  the  shock-boundary  layer  interaction)  is  generally  taken  to  be  less  than 

one.  Typical  values  of  Ax/Sq  for  shock-turbulent  boundary  layer  interactions 

are  0.17  in  Reference  16,  0.45  in  Reference  10  and  1.0  in  Reference  36. 

The  presence  of  flow  separation  has  an  important  influence  on  Ax/6  , 

requiring  it  to  be  less  than  0.5  in  general.  Experience  with  the  numerical 
1 8 

code  in  computing  inlet  flows  that  do  not  exhibit  flow  separation  has 

shown  that  the  criterion  Ax/6q  ;  1  is  generally  sufficient. 

The  coordinate  transformation  program  can  be  used  with  either  uniform 

or  arbitrary  mesh  spacing  in  the  x-direction.  For  uniform  mesh  spacing, 

the  procedure  is  to  estimate  the  values  of  <5  as  discussed  above,  using 

either  experimental  data  or  approximate  techniques,  and  to  choose  a  value 

of  Ax  that  insures  that  the  criterion  on  Ax/<5  is  met  at  all  streamwise 

o 

loca  tions . 

2 .  Boundary  Layer  Mesh  Spacing 

The  mesh  spacing  in  the  boundary  layer  or  n-direetion  is  influenced 
by  three  major  factors.  First,  there  must  be  sufficient  resolution  of  the 
boundary  layer  by  the  computational  sublayer  mesh  and  ordinary  mesh  combined. 
Typically,  a  total  of  10  to  15  ordinary  mesh  points  should  be  within  each 
boundary  layer  on  the  ramp  and  cowl  at  all  streamwise  locations.  Of  course, 
since  the  entire  flowfield  is  not  known  a  pri ori ,  this  achievement  of  this 
criterion  can  only  be  determined  a  posteriori.  The  practical  impact  of 
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this  criterion  for  the  purpose  of  preliminary  analysis  is  to  assist  in 

determining  the  appropriate  distribution  of  mesh  points  on  the  upstream 

and  downstream  boundaries,  using  program  BNDRY.  Past  results  in 

18 

calculation  of  high  speed  inlet  flows  have  shown  that  the  satisfaction 

of  the  requirement  of  10  to  15  ordinary  mesh  points  within  each  boundary 

layer  at  the  upstream  and  downstream  boundary  generally  insures  achievement 

of  the  requirement  at  all  streamwise  locations.  Insofar  as  the  computational 

sublayer  mesh  is  concerned,  a  sufficient  number  of  points  are  required  in 
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order  that  the  viscous  sublayer  be  adequately  resolved.  This  implies  ’ 
that  Ay„T  <  5  at  all  streamwise  locations,  where  A' 


=  AySL 

u  N/v  , 
*  w 

with 

r  /o 

and  v 

a 

p. 

it 

w  w 

V 

w  w 

r  ,  o 

w  w 

and  u 

w 

are  not 

known  a  priori,  and  in  practice  an  engineering  estimate  is  used  (see  below) 
to  insure  satisfaction  of  the  requirement  at  a  selected  number  of 
streamwise  locations. 

The  second  factor  influencing  mesh  spacing  in  the  boundary  layer 
direction  is  the  requirements  of  efficiency  and  accuracy  imposed  by  the 
use  of  the  computational  sublayer  technique.  As  discussed  in  Section  T.C.l, 
this  implies  (see  Figure  7) 


v  v 

20  —  <  v  ’  <  60  — 1 
u*  '  'm  '  u* 


(61 ) 


Assuming  3p/3n  -  0  in  the  boundary  layer,  the  length  scale  v  / u  is  given  by 


u  RT 
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p  "  U  .  , 

we  f  w 


j~  T-  (6  2) 

*  cr  1 
f  w 

where  c^-  =  t  /  [  (1/2)  (,>eU^)  ]  and  the  subscript  e  refers  to  conditions  at 


the  edge  of  the  houndarv  layer.  The  utilization  of  this  expression  requires 


estimates  for  (the  Mach  number  at  the  boundary  layer  edge)  and  the 

wall  pressure  p  ,  obtained  either  from  experimental  data  or  a  simplified 

inviscid  analysis.  The  value  of  T  (and  hence  p  )  can  then  be  estimated 

w  w 

(for  adiabatic  walls)  from 

T  =  (1  +  --  Pr  M2)T  (63) 

w  2  t  e  e 


where  Tg  is  approximated  by 

T  =  T  /[l  +  -(y  1)  M21  (64) 

e  t  z  e 

CO 

where  Tr  is  the  freestream  stagnation  temperature.  The  value  of  U 

'-oo  e 

may  be  obtained  from  U  =  M  /yRT  .  The  value  of  c,  mav  be  estimated  by 

e  e  e  f  7 
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a  variety  of  engineering  techniques.  It  is  important  to  note  that  cf 

can  be  increased  substantially  by  boundary  bleed.  Experience  in  computing 
1 8 

high  speed  inlets  indicates  that  typical  values  of  c^.  in  bleed  regions 
can  be  as  large  as  0.01  depending  on  the  bleed  mass  flux. 

In  practice,  it  is  generally  only  necessary  to  estimate  an  appropriate 
value  of  y^  from  (61)  at  a  few  points.  These  points  are  typically  taken 
to  be  upstream  and  downstream  of  each  shock -boundary  layer  interaction, 
since  the  value  of  v  / u .  varies  inversely  with  the  wall  pressure.  The 
above  expression  (62)  can  also  be  used  to  determine  the  number  of  points 
required  in  the  sublayer  in  order  to  satisfy  the  requirement  Ay^  <  5 
discussed  above. 

The  third  factor  influencing  mesh  spacing  in  the  boundary  layer 
direction  is  given  by  the  requirement  of  numerical  stability  in  the 
computational  sublayer  in  the  presence  of  bleed.  The  condition  is  that 
the  bleed  cell  Reynolds  number  Re1  defined  by 
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(65) 


Re-  = 
m 


I m!  Ay 


SL 


2(ti  +  e) 


be  everywhere  less  than  approximately  0.25  within  the  computational 
sublayer.  Noting  that  the  maximum  value  of  Re-  general Lv  occurs  on 

Til 


the  wall,  this  condition  can  be  approximated  by 


!ml  Ay 

— = - —  <  0.25  (66) 

2u 

w 

at  all  streamwise  locations.  Since  the  values  of  m  are  known,  and 
U  is  a  function  of  T  t  T.  for  adiabatic  walls,  this  condition  can 

W  W  too 

be  directly  evaluated  at  all  locations. 

In  some  instances,  there  is  a  boundary  layer  on  the  lower  surface- 
only  (e.g.,  mesh  A  in  Figure  8).  The  upper  boundary  n  -  1  must  therefore 
be  placed  at  a  sufficiently  large  distance  from  the  lower  surface  in  order 
that  the  freestream  boundary  conditions  are  applicable.  This  distance 
is  typically  taken  to  be  approximately  56,  where  6  is  the  boundary  layer 
thickness  at  the  downstream  end  of  the  mesh The  use  of  equations  (10) 
to  (32)  will  cause  a  decrease  in  mesh  spacing  near  the  upper  boundary. 

This  decrease  will,  not  serious!  v  affect  code  performance,  and.  lias  surprising! 

•k 

yielded  an  improved  solution  for  a  flat  plate  tuibulent  boundary  Liver. 


* 

In  particular,  calculations  for  a  fully  turbulent  flat  plate 
boundary  layer  have  been  performed  using  a  rectangular  computational 
domain  whose  spacing  in  y  is  symmetrical  about  the  line  midwav  between 
the  lower  and  upper  boundaries.  This  type  of  coordinate  system  can 
easily  he  obtained  by  choosing  JI,  to  be  odd,  DY0  -  1)Y1  and  rij  =  r;  -  0.5 
(  i  ,e  . ,  J00=J01=J10=J11=( JL+1 ) / 2 ;  set'  Sections  1  T.B  and  11  .(')  .  It  has 
been  found  that  with  this  type  of  coordinate  system,  there  is  no  char¬ 
acteristic  "hump"  in  6)  immediately  downstream  of  the  leading  edge. 

This  hump  had  been  observed  with  the  code  for  fuliv  turbulent  flat  plate 
boundary  layer  calculations  with  mesh  spacing  in  v  consisting  of  a 
uniform  stretching  within  the  boundary  layer  followed  bv  a  constant 
mesh  spacing. 
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In  Figure  9,  a  portion  of  the  flowfield  of  a  simulated  high  speed  inlet 

is  illustrated.  The  configuration  corresponds  to  the  upstream  inlet  portion 
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of  Case  35  of  the  MCAIR  study.  The  freestream  Mach  number  M  is  3.51, 

and  the  freestream  total  temperature  Tt  and  total  pressure  are  576°R 

and  7,099  lbf/ft2, respectively .  The  freestream  static  pressure  and 

2 

temperature  are  91.8  lbf/ft  and  166°R.  A  fully  developed  turbulent  bounder’ 
on  the  ramp  enters  a  converging  duct  formed  by  two  plates.  The  shock  wave 
formed  by  the  cowl  surface  deflection  intersects  the  ramp  at  x  =  20  in. 

The  values  of  the  boundary  layer  bleed  mass  flux  on  the  ramp  and  cowl  are 
indicated . 

First,  the  appropriate  streamwise  mesh  spacing  Ax  is  determined.  Due 

to  the  presence  of  substantial  boundary  layer  bleed  on  the  ramp  in  the 

vicinity  of  the  shock -boundary  layer  interaction,  it  is  unlikely  that  flow 

senaration  will  occur,  and  hence  the  criterion  Ax/ 5  '  1  is  satisfactory 

o 

on  the  ramp.  The  experimental  values  of  .5  on  the  ramp  vary  fro-  0.25  in. 
at  x  =  1.9.2  in.  to  0.20  in.  at  x  =  24  in.  as  indicated  in  Figure  9,  and 
thus  a  value  Ax  =  0.25  in.  is  sufficient.  Since  the  maximum  allowable 
value  of  IL  is  40,  this  places  the  downstream  boundary  at  x  =  23.63  in. 
upstream  of  the  reflected  shock  intersection  on  the  cowl.  Thus  the  cowl 
boundary  layer  is  not  subject  to  any  strong  viscous-invisri d  interaction 
within  this  computational  domain.  Based  on  the  experimental  value  of 
'  =  0.10  in.  at  x  *  25  in.,  this  yields  Ax/A  2.5  on  the  cowl,  which  is 
satisfactory  for  this  region. 

Secondly,  the  requirement  of  resolution  of  tin*  boundary  layers  is 
considered.  Based  on  the  experimental  values  of  the  boundary  Inver 


thickness  on  the  ramp  and  cowl,  the  following  criteria  are  established. 


Boundary  Surface  Criterion 


Upstream 

Ramp 

10 

to 

15 

points 

below 

y  =  0.25 

in . 

Cowl 

10 

to 

15 

points 

above 

y  =  2.09 

in . 

Downstream 

Ramp 

10 

to 

15 

points 

below 

y  =  0.20 

in . 

Cowl 

10 

to 

15 

points 

above 

y  =  0.88 

in . 

Ordinary  mesh  points. 


These  criteria  are  directly  utilized  in  determining  the  boundary  mesh 

distribution  using  program  BNDRY . 

Next,  the  computational  sublayer  requirement  (61)  is  considered. 

Considering  the  ramp,  it  is  evident  that  the  major  changes  to  the  boundary 

layer  arise  due  to  the  shock-boundary  layer  interaction  and  mass  bleed. 

It  is  sufficient,  then,  to  consider  the  ramp  at  two  locations,  e.g., 

x  =  13.88  in.  (upstream  boundary)  and  x  =  21 .8  in.  At  the  first  location, 

Te  =  166°R  and  =  534°R  from  (63).  Using  an  estimate  of  =  1.7  x  1.0  ^ 

(typical  for  M  =  3.5  adiabatic  flat  plate  boundary  layer  at  these  Reynolds 

numbers),  and  noting  that  pw  '  pn  ,  the  estimated  value  of  v  / uA  is 

3.3  x  10  ^  ft.  From  (61),  the  desired  range  of  v’  at  x  =  13.88  in.  is 

•  m 

-4  -  3 

approximately  6.5  x  10  to  2,0  x  10  ft.  At  the  second  location,  the 

experimental  Mach  number  M  is  2.5  (this  could  also  be  estimated  from  a 

simple  inviseid  analysis,  knowing  the  cowl  surface  deflection)  and  the 

wall  pressure  p  5  p  .  From  (63)  and  (64),  T  -  256°R  and  T  •  544°R, 
w  e  w 

Since  boundary  1  aver  bleed  is  present  at  x  =  21  .8  in.,  an  estimate  of 

8  x  10  is  used  for  c e  (based  on  experience),  yielding,  /u.,  =  4.4  x  10  ^  ft, 

indicating  a  desired  range  of  8.7  x  10  ^  ft  to  2.6  x  10  *  ft  fur  v ’  at 

m 

x  =  21.8  in. 


On  the  basis  of  this  information,  it  was  decided  to  use  y'  =  5 

m 

from  x  =  13.875  in.  to  approximately  x  *  19  in.,  decreasing  to  y*  = 

at  x  =  20  in.  and  constant  thereafter.  Although  the  upstream  value 

is  slightly  below  the  range  indicated,  it  nonetheless  does  not  i.rpa 

code  efficiency  because  of  the  smaller  values  of  v'  for  x  >  30  in. 

'  m 

values  of  v ,+  are  thus  expected  to  varv  from  15  at  x  =  13.83  in.  to 
'  m 

57  at  x  =  23.63  in. 

The  number  of  mesh  points  in  the  computational  sublay. -r  on  tin- 
determined  by  the  requirement  that  Ay*  f  5  everywhere,  fin  -he 
solution  is  extremely  rapid,  the  maximum  allowable  number  of  oo'nis 
20)  may  be  employed.  For  20  points  in  the  sublayer,  it  is  evident 
Ay*  varies  from  0.79  to  3.0  on  the  ramp,  which  is  enti  rol  v  sa*  i  sr  : 

oLi 

Finally,  it  is  necessary  to  check  that  the  condi  tier.  ffti.)  is 
everywhere  in  the  ramp  boundary  layer.  In  bleed  zone  1,  m  -  -'<.85 

2  -  -i 

si  ugs/f  t^-sec  and  the  maximum  value  of  -\y  is  3.0  10  f  t.  /  !  9  "r 

2.6  x  10  5  ft,  yielding  a  value  ;ni  Av  /2  0.1  1.  In  Meed  zone 

'  Mj  V 

—  s 

the  maximum  value  of  Av  is  1.3  x  10  ft,  and  .  t. .  ‘v  /.'  ,  0.11 

'SI.  SL  V 

values  are  within  the  limit  indicated  by  (66). 

An  analysis  similar  to  the  above  Is  also  required  'or  the  owl 

results  indicate  that  v.  „  -  v.  „  ,  should;  behave  similar  to 

'  i,  JL  -  i  ,  JT.-l 
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B.  Determination  of  Mesh  Spacing  on  Upstream  and  Downstream  Boundaries 
Using  Program  BNDRY 

1 .  Introduc  tion 

The  program  BNDRY  is  employe'  to  determine  the  distribution  of 
mesh  points  in  general  on  the  upstream  and  downstream  boundaries.  These 
boundaries  are  assumed  to  be  vertical  lines,  with  a  mesh  point  distribution 
given  by  equations  (30)  to  (32).  The  purpose  of  the  program  is  to  determine' 
the  values  of  the  parameters  C^,  C^,  and  n0  that  yield  the  desired  mesh 
distribution  on  the  boundaries.  This  information  is  then  utilized  by  the 
coordinate  generation  program  in  determining  the  entire  mesh. 

As  discussed  in  Section  III. A,  the  following  information  is  required 
for  each  boundary:  (1)  the  desired  height  of  the  first  ordinary  mesh  point 
above  the  lower  and  upper  boundaries,  i.e.,  y,  _  -  v.  and  v.  -  v.  , 

1)Z  X  y  i.  1  y  L<  I  ,  ij  l.  .1 

respectively  (where  i  =  1  or  IL) ,  (2)  the  boundary  layer  thickness  on  the 

lower  and  upper  boundaries,  (3)  the  height  of  the  boundary  y,  -  v.  . 

1  y  J  L  1,1 

All  dimensions  in  the  program  are  in  feet.  The  user  chooses  a  range  of 

values  of  and  to  be  investigated.  For  each  value  of  n  and  q0,  the 

program  computes  by  Newton's  method  the  values  of  C.^  and  C?  from  (34)  and 

(35)  that  yield  the  specified  values  of  v.  _  -  y.  ,  and  v.  ,,  -  v.  _  ., 

■  j  ,  2  '  i ,  1  ■  l. ,  JL  '  i ,  JL-  i 

and  prints  the  entire  mesh  distribution  on  the  boundary  y  .  ,  j  =  1,  .... 

1  t  J 

JL  (where  i  =  1  or  IL) .  The  user  then  chooses  the  particular  set  of  values 
for  and  that  yield  the  desired  number  of  ordinnrv  mesh  points  in  each 
boundary  layer  (c.g.,  10  to  15). 

The  following  table  provides  information  on  the  resources  require:' 
to  execute  program  BNDRY  on  the  OYBLR  175. 
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TABLE  1.  RESOURCES  REQUIRED  FOR  EXECUTION  OF  PROGRAM  BN' DRY 


Resource 


Details 


Computer  time  <  5  seconds  (typical) 

Input/Output  time  <  2  seconds  (typical) 

Core  Memory  <  60,000  words  (octal) 

Files  required  INPUT,  OUTPUT 


2 .  Description  of  Input  Variables 

The  input  to  program  BNDRY  consists  of  two  cards  (card  "ima 
or  Lines  in  the  case  of  batch  input  from  a  terminal).  'Hie  Juf  nition 
format  of  the  input  data  is  indicated  in  Table  2.  For  all  source  •:  ed 
variable  names,  the  symbol  0  implies  zero,  and  the  svmbol  0  implies  ' 
letter  0. 

If  the  estimates  for  C,  and  are  t far  'if,  the  program  mar- 
verge  to  CjL  —  C'2  ~  0  (and  noss  ibl  v  display  an  error  message  indi.  satin 
underflow  and/or  overflow  has  c,.  .rtvdl  .  s  ince  0  anci  C0  must  .-•<  pi¬ 
th  e  program  should  then  be  rerun  v  i  th  a  i  t<-rnate  eho i cos  for  and  c 
With  some  experience,  this  difficult'.-  will  be  avoided. 


TABLE  2.  INPUT  DATA  FOR  PROGRAM  BNDRY 


Line  1: 

DY0,  DYI,  RH , Y0 ,  CEST1, 

LEST  2 

(Format  6F10.4) 

Fortran 

Definition 

Pumge  or  Value 

DY0 

yi,2  -  yi,l 

Units:  feet 

DYI 

y  i,JL  i, JL-1 

Units:  feet 

RH 

yi,JL  "  Yi,l 

Units:  feet 

Y0 

yi,l 

Units:  feet 

CEST1 

Rough  estimate  for  Cl 

Tv pie a  1 1 v  b etwe e n 
1  and  U) 

CEST2 

Rough  estimate  for  C2 

ly pi  rally  between 
1  and  10 

Line  2: 

JL,  J0STAR,  J0END,  J1STAR 

,  JIEND  | 

(Format  5T5) 

Fortran 

Def ini tion 

Range  or  Value 

JL 

Number  of  points  in 
tj— di  rect  ion 

J0STAR, 

Define  smallest  and 

2-  J0S1AK 

J0END 

largest  values  of  >] ^ 
to  be  investigated. 

J0STAR-  J0END 
J0END-  JL-2 

i.e.,  (  J0STAR-1.)  / 
(JL-1)  <  1'; j 
(J0KND-1) /tJL-1) 

.]  1STAR, 

Same  as  above  bur  for 

S  i •  '  a  r  t  o  ah  ove 

J1END 

3. 


Source  Code  Notation 


The  definition  of  the  pertinent  source  code  variables  is  given 


below.  For  any  array,  the 

index  J  is  equivalent  to 

the 

subscript  j  (erg 

ETA(J)  is  equivalent  to  ru) 

• 

Fortran 

Definition 

Range  or  Value 

DETA 

An  =  1/ ( jl-1) 

DY 

»<•»,-  n.s-  n.1-1 

(l  =  1  or  1L) 

l-nits:  feet 

ETA1 

ni 

ETA2 

n2 

YY 

YY(J)  =  y  .  (i  =  1  or 
i>  1 

1 L) 

Anils:  foe*: 

4.  Sample  of  Output:  Upstream  Inlet  Region  for  MU 

\i_R 

Case  >5 

A  portion  of  the  sample  output  of  prog-. rat, 

AY  for  the  crier: 

of  the  upstream  inlet  region  of  the  MCA1R  inlet  dim 

efi  ?  n  f  ee  t  i  on  •  -  i 

is  indicated  in  Figure  10. 

This  calculation  was  ;><.  - 

■  f  i  r 

.re)  as  :  I>f4  r‘U' 

parameters  presented  in  the 

example  in  the  previous 

S  O  i. 

tl  o,r  As 

earlier,  the  requirements  for  the  upstream  hnun.'.-rv 

.>f 

•A.  ■-■V.V!  ■  ■ 

region  are  the  following: 

Requirement 

Cr 

i  t.  e  r  i  n  n 

Resolution  of  ramp  boundary 

layer  a.  .10  :■ 

4  '  S 

1  )  s 

anil  nary  no  Mr  -  ' 

h  ’  v  !  .  ’ 

.  .  :».o  >;  j * 

i.  -* ; 

Resolution  of  cowl  boundary 

layer  a.  10  t . 

o.vrnirv  no 

h .  v 

■  L . .  1 1 

- 

V  t  .  •  \I)  V,  ' 

It  is  evident  from  Figure  10  tluiL  values  of  •• .  (  )  ,'t  0  - !  >  | 

and  i|  =  (Jl-1)  /(JL-1)  =  27/47  are  sat  is  f  ac  torv  .  liu  e.  i  i  ’lid  i  v.: 
of  C  and  C?  are  2.98068  and  2.87725,  respectively.  Tin  1  ■  ui  p.ir.mvt  e? 
.11,  C  ,  C  are  utilized  in  the  coordinate  general  ion  ;>r  n  ran. 
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C.  Generation  of  Coordinate  Trans  format  i  on  Using  Program  too  kb 
1 .  Introduction 

The  program  COORD  is  employed  to  calculate  the  i oord : note  tran 
formation.  Prior  to  the  execution  of  this  program  the-  desired  streamwise  and 
boundary  layer  mesh  spacing  must  be  determined  by  a  preliminar;  analysis  as 
discussed  in  Section  III. A. 

In  the  present  form,  the  COORD  program  incorporates  the  general  category 

33 

of  inlet  shapes  employed  by  the  MCAIR  study,  including,  o 1  course,  the  flat 
plate  portion  of  the  ramp  upstream  of  the  inlet  entrance,  using,  uniform  strea 
wise  mesh  spacing.  The  program  can  be  simply  modified  lo  handle  arbitrary  in 
let  geometries,  as  detailed  in  subroutine  IN IT  of  the  program. 

The  following  table  provides  information  on  the  resources  required  to 
execute  program  COORD  on  the  CYBER  175.  Care  must  be  taken  to  insure  that 
sufficient  computer  time  is  requested  in  order  to  complete  the  number  of 
iterations  specified.  The  program  permits  restarting  and  continuation  of  the 
calculation  if  convergence  is  not  achieved  within  the  requested  number  of 
iterations . 


TABLE  3.  RESOURCES  REQUIRED  KOR  EXECUTION 
OF  PROGRAM  COORD 


Resource 

Computer  time 
Input/Output  time 
Core  Memory 
Fi  l.es  Fequ  i  red 


De La  i  is 

•bO  seconds  (tvpical) 
■  10  seconds  (,  typical.) 
150,000  words  (octal) 
INPUT,  nriTLT,  XY, 
OVR1.AP,  KS’iART 
(See  See  l '  -'ll  11  T  .  (' .  5) 


2 .  Description  of  Input  Variables 

The  input  to  program  COORD  consists  of  a  variable  number  of  cards. 


The  definition  and  format  of  the  input  data  is  indicated  in  the  required 


order  in  Table  4.  In  those  cases  where  the  definition  of  an  input  variable 
is  lengthy,  the  reader  is  referred  to  Section  III.C.4.  Line  numbers  are  not 
assigned  since  the  total  number  of  input  lines  is  variable.  Additional  in¬ 
formation  is  given  in  Section  III.C.4.  Useful  guidelines  for  the  selection 
of  certain  input  parameters  are  provided  in  Table  5. 
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TABLE  4.  INPUT  DATA  FOR  PROGRAM  COORD 


Line: 

M,  N,  MAXITR,  IRDZET,  IOVLP, 

I  RETRY 

(Format  615) 

Fortran 

Definition 

Range  or  Value 

M 

Number  of  points  in 
C-direction  (i.e.,  IL) 

<40 

N 

Number  of  points  in 

<48 

H-direction  (i.e.,  JL) 

MAXITR 

Maximum  number  of 

_>1. 

iterations  of  equation 
(24) 

IRDZET 

See  Section  III. 0.4 

0  or  ] 

IOVLP 

See  Section  IJI.C.4 

i.<:rovLP<  m 

I RETRY 

See  Section  III.C.4 

0  or  1 

- 1 


Line: 

J00,  J01,  J10,  Jll, 

IRE  11 , 

IRLF2  | 

(Format  615) 

Fortran 

Definition 

Knnr.o  or  V 

1  1  U 

J00 

Defines  n  at 

C  =  0 

1 '  100  -  .  1  1  0 

J01 

Defines  rt _  at 

=  ] 

2-M01-J  11 

J10 

Defines  n.,  at 

4  =  0 

.  00  .11 0--.IL 

Jll 

Def  Lnes  rjn  at 

-  ] 

.1  I'/'  *5  i  !  . 1  i . 

1REF1 

See  Section  II 

I  .  C  .  4 

1-  I  i  I  -  M 

I  REF  2 

See  Sect i on  I I 1 . C . 4 

1-;. 
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TABLE  4.  CONT'D 


Line: 


(Format  215) 
Fortran 
JREF 
IDSCNT 


Line: 


(Format  1615) 


Line: 


(Format  415) 
Fortran 
ID(K) , 

K=l,  . ..,  20 


Line: 


(Format  1615) 


Line: 


(Format  415) 
Fortran 
JD(  K) 

K=l,  ....  20 


JREF,  IDSCNT 


Definition 

See  Section  III.C.4 

See  Section  III.C.4 


Range  or  Value 

1<JREF<N 

0_<IDSCNT<20 

A 


ID(1) ,  ID(2) ,  ....  ID( 16) 


ID(17) ,  ....  ID(20) 


Definition 


See  Section  III.C.4 


V  es  required 


J 

Range  or  Value 
KID(K)<M 
K=1 ,  ....  20 


JD(1),  JD( 2) ,  ...,  JD(16) 


JD(17) ,  ...,  JD(20) 


Def inition 


See  Section  III.C.4 


as  required 


J 

Range  or  Value 
JD(K)  =  1  or  S 
K=l,  ...,  20 


•  If  IDSCNT  =  0,  these  lines  are  omitted 


TABLE  4.  CONT'D 


Line:  DY00,  DY01 ,  DY10,  DY 1 1 

(Format  4F10.4) 


Fortran 


Def i ni t i on 


Ranee  or  Value 


DY00 

DY01 

DY10 

DY11 


V1 >  2  -  yi,l 

Units: 

feet 

Y IL  ,  2  ~  y  IL  ,  1 

I'n  i  ts : 

feet 

yl,JL  ~  Y1 , JL-  1 

1'n.  i  ts : 

I  eet 

r- ( 

>- 

1 

hJ 
1— 1 
>b 

Un its: 

f  eel 

Line:  CC1(1),  CC2(1),  CC1(M), 


CC2( M) 


(Format  4F10.4) 

Fortran 

CC1(1) 

CC2( 1) 

CC1(M) 

CC2(M) 


Def i nit  ion 
Value  of  Cl  at  (  =  0 

Value  of  C2  at  r,  -  0 

Value  of  Cl  at  ",  =  I 


R_an_ge  or  Value 


\ 


F  rum 

l'rograr 

UNDRY 


Value  of  C2  at  -  1 


Line:  SCALEl,  SCALE 2,  SCALE.  1,  SCAi.F.4 


(Format  4F10.4) 

Fortran  Def  in i  t  i on  Rar.jtv  or  VaJ_ue 

SCALE1  Provides  smooth  o.i/4 

transition  for 
T(4,n)  between 
three  regions  in 
equation  (2). 

SCALE2  Same  as  above  0.174 

SCALE3  See  Section  Ill.C.-* 


SCALE4 


See  Section  II T . C  .  4 
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TABLE  4 .  CONT ' D 


Line: 

S0(1),  Sl(l),  CONVER 

(Format  3F10.4) 

Fortran 

Definition 

Range  or  Value 

S0(  1) 

See  Section  III.C.4 

Units:  feet 

Sl(l) 

See  Section  III.C.4 

Units:  feet 

CONVER 

Convergence 

Criteria  (See  (28)) 

See  Table  5 

**  Line: 

X0,  DX,  XH1,  DELTAC 

(Format:  4F10.4) 

Fortran 

Definition 

Range  or  Value 

X0 

Value  of  X  at 
upstream  boundary 
of  mesh  region 

Units:  feet 

DX 

Mesh  spacing  in 
X-direction 

Units:  feet 

XH1 

Height  of  inlet 

throat 

Units:  feet 

DELTAC 

Cowl  angle 

Units:  Degrees 

**This  applies  to 

MCAIR  inlet  (Reference 

33) 

Note  that  X  =  0 

at  the  ramp  leading  edg 

e  f or 

the  MCAIR  inlet. 
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TABLE  5.  GUIDELINES  FOR  INPUT  DATA  FOR  PROGRAM  GOURD 


Parameter(s) 

J00,  J01 

J10,  Jll 
IREF1 

IREF2 

JREF 

IDSCNT,  ID(K), 

DY00,  DY10 

DY01,  DY11 
CC1(1),  CC2(1) 

S0(1) ,  SI ( 1) 

CONVER 


Guideline 


Use  same  value  fur  both,  in  general,  for 
all  overlapping  mesh  regions  (e.g.,  meshes 
B  and  C  in  Fig.  8) 

Same  as  above 

Care  must  be  taken  that  si  pep  corresponds 
to  same  phys i cal  location  when  Overlapping 
mesh  regions 

Same  as  above 

Same  value  required  for  all  overlapping 
mesh  regions 

JD(K)  Care  must  be  taken  to  insure  that  e.e  n 

surface  slope  discontinuity  is  a  j  wavs  referred 
to  its  proper  physical  location  wh>n  overlapping 
mesh  regions 

When  overlapping  mesh  regions,  the  values 
of  DY00  and  DY11  can  fit  obtained  frort  the 
coordinates  of  the  upstream  mesh  ti.o 

restart  station 

These  values  are  always  obtained  1  rota 
program  BNDRY 

When  overlapping  mesh  regions,  uie.»c  values 
can  be  obtained  from  the  output:  of  (he 
upstream  mesh  at  the  restart  station 

These  values  represent  the  arc  dist.ar.ee 
along  the  lower  and  upper  stir  faces  from  the 
leading  edges  on  the  respective  surfaces, 
even  though  the  leading  edge  does  not  lie 
in  that  particular  mesh  region 

Usual  values  are  If!  K  or  less.  Typical  ly, 

CONVER  is  100  to  1000  times  smaller  than 
the  smallest  mesh  spacing  in  r  eel  in  either 
r.  or  M  direction. 


3.  Flow  Chart 


The  flow  chart  for  program  COORD  is  indicated  below.  The  names  of 
the  subroutines  called  are  shown  by  capital  letters  in  parentheses. 


Subroutine 

INIT 


BECIN 


CONTINUE 


1 
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4 .  Source  Code  Notation 

The  definition  of  the  pertinent  source  code  variables  is  given  below. 
For  any  array,  the  index  I  is  equivalent  to  the  subscript  i  (e.g.,  the 
Fortran  variable  ZETA(I)  is  equivalent  to  4.  in  that  ZETA(l)  =  4^, 

ZETA(2)  =  42.  etc.)  and  the  index  J  refers  to  the  subscript  j. 


Fortran 

CC1(I) 

CC2(I) 

CONVER 

DELTAC 


Definition 

: k 

Value  of  C^  at  4.  (see  (2)  and  Note  1  ) 
Value  of  C2  at  4.  (see  (2)  and  Note  1) 
Convergence  criteria  (see  (28)) 


Cowl  angle  for  MCAIR  inlet  (in  degrees) 


DX 

DY00 

DY01 

DY10 

DY11 

EDA(J) 

H 

HH 

ID(K) , JD(K) 


IDSCNT 


Mesh  spacing  in  x-direction  (in  feet) 

Value  °f  Yj  ■)  '  Y\  1  (in  feet) 

Value  of  yIL>2  -  yILJ  (in  feet) 

Vaiue  of  y1?JL  -  y^^  (in  feet) 

Vaiue  of  y^^  -  (in  feet) 

nj 

A4 

(AC)2 

Denote  locations  on  the  lower  and  upper  boundaries 
where  surface  slope  is  discontinuous.  Denoting 
ID(K)  by  l  and  JD(K)  by  m,  the  K*^1  discontinuity 
lies  between  r,  f  and  -,r+|  on  the  '  i-ie  pm.  A 
maximum  of  20  discontinuities  is  allowed.  The 
discontinuities  may  be  listed  in  any  order.  See 
Note  2. 

Total  number  of  discontinuities  in  surface  slope 
on  lower  and  upper  boundaries  (maximum:  20). 
Corners  are  not  included.  If  there  are  no 
discontinuities,  set  IDSCNT  =  0. 
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IOVLP 


1 


Used  in  generating  overlapping  r:usti  regions, 
such  as  regions  li  and  C  in  Figure  8.  The  curve 
5  =  CTnvl  p  °f  the  upstream  mesh  is  taken  to 
represent  the  t,  =  0  boundary  of  the  downstream 

mesh . 


IREF1 


IREF2 


IRDZET 


I RETRY 


ITER 


In  conjunction  with  the  quantity  SCALE3,  this 
term  defines  the  streamwise  extent  over  which 
the  height  of  the  first  row  of  ordinary  mesh 
points  above  the  lower  surface  (i.e.,  2  "  Yi  ]) 

changes  from  DY00  to  DY01  .  See  Mote  3. 

In  conjunction  with  the  quantity  SCALF.4,  this  term 
defines  the  streamwise  extent  over  which  the 
distance  of  the  first  row  of  ordinary  mesh  points 
below  the  upper  surface  (i.e.,  vj  ji  -  y ^ 
changes  from  DY10  to  DY]  1  .  See  Note  3. 

Used  in  generating  overlapping  mesh  regions.  If 
IRDZET  equais  1,  the  upstream  mesh  is  read  in 
order  to  define  the  -  0  boundary  of  the  downstream 
mesh.  If  mesh  overlapping  is  not  to  be  performed, 
set  IRDZET  equal  to  0. 

Used  to  continue  iterating  a  mesh  that  did  not 
converge  in  the  previously  specified  number  of 
iterations.  If  I RETRY  equals  1,  the  last 
results  are  read  and  the  iteration  continued. 

If  this  restart  feature  is  not  desired,  set 
IRETRY  equal  to  0. 

Denotes  iteration  number  during  the  course  of 
the  calculation. 


JREF  Used  in  defining  array  XNU..I)  (see  below). 

For  J  JREF.  the  quant  i  tv  XN(T,J)  is  the  normal 
distance  of  the  point  t  r. ;  j ,  y  j  j)  above  q  =  0. 
For  J  >  JREF,  tin  quant i t  v  XN(1,J)  is  the  normal 

distance  of  (x.  .,v.  .)  from  q  =  1  . 

1,1  i.l 

Note:  In  utilizing  program  INI. FT.  it  is 

imperative  that  the  id  levins  he  observed: 

IP!  !  0  J!M:F 
JR!  FI  JIF  F 


where  JREK0  mJ  IF!  ’  1 :  • 

program  INl.fF. 

ra  r  i  ah  i  t's  defined  in 

J00 

Defines  value  el  ■  1  at 

n  =  (J00-D/C  11. -n 

11  (see  (.)))  through 

J0 1 

Defines  q^  at  -  ! 

.110 

Defines  at  ;  ' 1  tiiro1 

igh  ,  -  ( JI0-1  )  /(.n.-i ) 

J1  1 

Defines  ■  ,  at  ! 
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K  An 

KK  (An)2 

M  Number  of  points  in  5 -direction  (i.e.,  IL)  . 

Maximum:  40 


Ml 


M  -  1 


MAXITR 


Maximum  number  of  iterations  of  algorithm  (24). 
See  TRETRY. 


N 


Number  of  points  in  n-direction  (i.e.,  JL) . 
Maximum:  48 


N1 

S0(I) 


Sl(l) 

SC AL El, SCALE 2 


N  -  1 

Distance  of  i*"*'  point  on  n  =  0  from  leading 
edge  of  lower  surface  (in  feet).  Negative  values 
imply  points  ahead  of  leading  edge.  The  distance 
is  measured  along  the  n  =  0  curve  (not  necessarily 
the  x-direction)  .  Note  that  the  value  cf  S0(1) 
is  positive  if  the.  4=0  boundary  of  a  particular 
mesh  region  lies  downstream  of  the  leading  edge 
of  the  lower  surface. 

Same  as  above,  except  refers  to  the  upper  surface. 

Provides  smooth  transition  for  T(r,n)  between 
the  t.iree  regions  in  equation  (2).  Typical 
value  is  0.374  for  each. 


SCALE3,  SCALE4  S-ee  IREF1 ,  TREE  2  and  Note  2 


T(I, J) 

T  (  C  i  , 

n j )  (see  (2)) 

TRAN ( I , J , 1 ) 

5*ifj 

(ZETAX(T, J) ; 

units:  it 

TRAN ( I , J , 2 ) 

(ZETAY(I,J); 

r  -U 
un its:  ft  ) 

TRAN  ( I , J , 3) 

(ETAX(1,J)  ; 

r  “K 

units:  ft  ) 

TRAN(T , J,4) 

nYi  .  « 

( F.TAY  ( I ,  J  )  ; 

r  -U 

un its:  ft  ) 

VARY 

Maximum  value  of  r 

(see  (27)) 

i  ,  J 

VARZ 

Maximum  value  of  r 

v-  •  (see  (27)) 

x  .  .  (inf ec  t ) 

1  ,  1 


V) 


X(I,J) 


XHl 

Height  of  throat  in  MCA  1  it  A 

XN(1,J) 

See  JREF .  Units:  f.-et 

X0 

Value  of  x  at  upstream  hoini' 
(MCAIR  inlet)  (in  feet),  Ti 
measured  from  the  ramp  lead 

Y(I,J) 

v.  .  (in  feet) 
i,  J 

ZETA(I) 

0 

Note  1:  The  values  of  CC1(1),  CC1(M),  CC2(1)  ami  <r2  0-t»  nv  r.tuird 
date  input.  They  are  found  from  program  P.N''!;'-  ir.-.i  from  ou  Up 
from  upstream  mesh  in  the  ease  of  overlapped  >1.  ’>“".ons. 

Note  2:  As  an  example,  consider  the  following  e:.nt  ..ss  N.n  tor  the  tow 
boundary : 

0  £  x  £  1:  y  ^  0 

1  <  x  £  2:  y=(x-l)sin  )0" 

2<_x<_3:  y  =  0  .  S 

Assume  mesh  points  were  distributed  on  I  lie  ■  -u:  ‘  «.-«•  at 

x  =  (n  -  1)  Ax,  n  =  1,2,  •••,  where  ‘  x  -  i.  ;  h,  ,  the 

following  values  would  represent  the  di.vont  in-  A  -  o 
x  =  1  and  x  =  2: 

ID(1)  =  4,  JD(1)  =  1 

ID( 2)  =  7,  JD(2)  =  1 


Note  i: 


The  height  of  the  first  row  of 
lower  surface  (i.e..  v.  »  -  v 


yi ,  2 


yi,2 


"  yi,l 


i  ,  1 

=  0 . 5(I)Y00  +  HY0  I  ) 


irdinarv  me 
)  is  S 1  a  i  j 


4-  O.5CDY0!  -  DY(/i0  )  t  a  nh  |  ( 

Thus,  the  height  changes  from  i)Y00  to  1  )s'0 ! 
range 


-2.65  Ar,*SCALE3  <  A  -  '  >  •••<.' 

For  example,  a  typical  value  <>l  Si 'At.?  t  i  .  .  >  , 

the  height  to  change  from  1»Y00  te-  l'Y0  i  >  . -r  . 

The  physical  extent  over  which  the  v.u  i  1 1  i .  ■  u  oi 

by  the  value  of  TRF.Fl  and  the  known  i  e  I  a  '  u 

and  r, .  on  the  lower  boundary, 
l 

A  similar  expression  holds  (or  the  uppoi  ui: ■ 


yi,.n,  Vi,JI.-i 


O.r)(l)Y10  4  PYI  1  1 
4-  0 . 5  ( DY  I  i  !  '■  !  0  I  • 


*.•  jpprnxi 


}  i*fl  iMUSt* 


}  l>  i  > 


et;  ion 
!  s 


for 

At 

.‘V 


*  Si  A], 
nu  !  - 


t  s  . 

riv  i  no- 
v  • ) 
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File  Structure 


OVRLAP 


RSTART 


Input  file  for  coordinate 
transformation  for  upstream 
region.  Used  in  generating 
overlapping  meshes 

Input  file  for  non-converged 
results 


Data  Structure 

a)  Convergence  achieved: 

X,  Y,  ZETAX,  ZETAY,  ETAX, 
ETAY,  XX,  S0,  SI 

b)  Convergence  not  achieved: 
X,  Y 

Same  as  a)  above 


X,  Y 


The  data  are  written  unformatted  in  the  sequence  indicated,  with  ail 
variables  except  S0  and  SI  having  array  size  40  by  48.  The  arrays  S0 
and  SI  are  each  of  size  40.  The  data  sequence  for  file  XY  is  thus 
X(l,l),  X(2 ,1) ,  X(40,l) ,  X( 1 , 2) ,  X( 2 , 2) ,  •••,  X(40,2),  •••, 

Y(l,l),  Y(l,2),  •••,  etc.  The  utilization  of  the  files  is  as  follows: 


File 

Type 

Utilized 

XY 

Output 

Required 

for  every  code 

exeout it 

OVRLAP 

Input 

Required 

only  if  IRDZET 

=  1 

RSTART 

Input 

Requi red 

only  if  I RETRY 

=  1 

6 .  Sample  Calculation:  Upstream  Inlet  Region  for  MCA!  R  Case  _  i_4 

In  Figure  11a,  the  coordinate  system  for  the  upstream  inlet  region  of 
MCAIR  Case  35  is  shown  (only  odd-numbered  rj-lines  a  tv  shown).  The  configura¬ 
tion  is  also  indicated  in  Figure  9.  The  input  data  is  indicated  in  Table  6 
and  the  input  variables  are  explained  in  Table  4  anti  Section  lil.C.4. 

It  is  Instructive  to  note  how  the  preliminary  anal vs  is  ft  the  upstream, 
inlet  region  is  used  in  the  input  to  program  COORb,  recai 1  lug  the  disc  ass  ion 
of  Section  III. A.  The  Wo  major  areas  of  the  pre.  1  ir  inary  ana  i s  i  s  are  dis¬ 
cussed  below: 

a)  Streamwlse  mesh  spacing 

The  chosen  value  of  Ax  =  0.25  in.  =  .0208  ft,  im!  a: -pears 
as  DX  in  line  10. 

b)  Boundary  Layer  mesh  spacing 

The  results  of  the  analysis  and  program  B'.'D'.Y  ini:  2 

J1  =  28,  Cl  =  2.89068  and  C2  =  2.87725  m,  the  upst  i  ■  a  -  h  . 

These  values  are  input  as  J00,  .110,  CC1  (1)  and  CC.  ,  n-sp.-,  ;  :  v.  1  . 

A  second  analysis  (not  detailed  in  this  report  i  i.i.ied  aid  -  . 

Jll  =  28,  CC1(M)  =  2.75649  and  CC2(M)  -  2. 7a  3/0.  du  st'  va.ues 
appear  on  lines  2  and  7. 

The  preliminary  analysis  also  indicates  th  u  v 1  =  5  IQ  ft 
at  the  upstream  boundary  on  the  ramp  and  cowl  ;)Y00  =  I’YK' 

.0005)  and  v'  =  2.5  x  10  ft  at  the  downs t  s  e.r  •  i>-.  und..rv  on  the  ram;- 

'  m 

and  cowl  (  i.e. ,  DY10  =  I)Y  1  1  =  .00024).  Tit  is  act  a  a  an  1  in.-  6. 

As  discussed  in  note  3  nf  Section  III.'  .  .  value  1.41  t  or 

SCALE  3  and  SCALK4  is  typical  (line  8).  A.-  in-ii  •  it .  a  : ::  Section 

III. A. 3,  it  is  desired  for  y  ’  to  decrease  t  r.v:  >.i:  \  !  ■’  t-  2.4  |'i 

m 

ft  in  the  vicinity  of  X  =  19.5  inch  on  hot  n  the  i  i  p  .nd  • ow 1  . 
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TABLE  6.  INPUT  DATA  FOR  COORDINATE  SYSTEM  OF 
UPSTREAM  INLET  REGION  (MCAIR  CASE  35) 


[Line  1 | 

M 

N 

MAXITR  IRDZET 

IOVLP 

(Format 

:  615) 

40 

48 

1000 

0 

0 

Line  2 

] 

m 

J$1 

J10 

Jll 

IREF1 

(Format 

615) 

18 

18 

28 

28 

23 

Line  3 

JREF 

IDSCNT 

(Format 

215) 

28 

1 

f 

Line  4 

1 

ID(1) 

( Fo  rma  t 

1615) 

5 

Line  5 

JD(1) 

(Format 

1615) 

48 

Line  6 

DY00 

DY01 

DY 10 

(Format 

4F10.4) 

.0005 

.00025 

.0005 

Line  7 

] 

CC1(1) 

CC2(1) 

CC1(M) 

(Format 

4F10.4) 

2.89068  2.87725 

2.75649 

Line  8 

SCALE1 

SCALE 2 

SCALE  3 

(Format 

4F10.4) 

0.374 

0.374 

1.5 

Line  9 

S0(1) 

Sl(l) 

CONVER 

(Format 

3F10.4) 

1.156 

-.093 

.000001 

Line  10 

] 

X0 

DX 

Xill 

(Format 

4F10.4) 

1.156 

.020833 

.06650 

I RETRY 
0 

IRF.F2 

23 


DY11 

.00025 

CC2(M) 
2. 74370 

SCALE4 

1.5 


DELTAC 

7.79 
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Noting  that  the  upstream  boundary  is  at  x  =  13. 88  inch  (Figure  9) 
and  the  streamwise  mesh  spacing  is  0.23  inch,  we  conclude  that  the 
station  x  =  19.5  inch  is  midway  between  the  stations  )  -  23  and  I  =  24. 
Thus,  we  take  I RE FI  =  IREF2  =  23  (line  2). 

In  the  following,  each  of  the  data  input  lines  in  Table  6  is  described. 

a)  Line  1 

A  mesh  of  40  points  by  48  points  (M  -  40,  N  =  48)  is  chosen. 

The  maximum  number  of  iterations  permitted  is  MAX1TR  =  1000.  The  up¬ 
stream  boundary  is  at  the  inlet  entrance,  and  thus  there  is  no  need 
to  read  a  prior  coordinate  system  (1RDZET  -  0,  I.0VL1’  =  0).  This  job 
represents  the  initial  attempt  to  solve  equations  (22),  and  thus 
IRF.TRY  =  0  (the  equations  (24)  will  converge  in  this  case  is  less 
than  1000  iterations) . 

b)  Line  2 

The  choice  of  the  values  for  the  variables  in  this  line  wore 
discussed  previously. 

c)  Line  3 

The  value  ofJREF  is  chosen  to  insure  that  (i  =  1,  ...,  M) 

always  lies  between  the  ramp  and  cowl  boundary  layers.  Based  on 
Figures  9  and  10,  we  choose  JREF  =  28.  We  also  note  from  Figure  11a 
that  there  is  only  one  discontinuity  on  the  lower  and  upper  boundary, 
namely  at  the  cowl  leading  edge;  thus  1DSCNT  -  I . 

d)  Lines  4  and  5 

From  Figure  9  and  the  chosen  value  of  Ax  =  0.24  inch,  we  note 
that  the  cowl  leading  edge  discontinuity  lies  on  the  upper  surface 
(j  =  N  =  48)  between  the  fifth  and  sixth  vertical  coordinate  lint' 

(i.e.,  between  1  =  5  and  i  =  6) .  Therefore,  we  choose  10(1)  =  8  and 
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JD( 1)  =  48. 


e)  Lines  6,  7  and  8 

The  choice  of  the  values  for  the  variables  in  these  lines 
were  discussed  previously. 

f)  Line  9 

As  indicated  in  Figure  9,  the  mesh  point  (i,j)  =  (1,1)  is 
located  at  13.88  inch  (1.156  it)  from  ramp  leading  edge;  thus 
S0(1)  =  1.156.  Similarly,  the  cowl  leading  edge  is  1.12  inch 
(.093  ft)  downstream  of  the  mesh  point  (i,j)  =(1,48)  and  thus 
Sl(l)  =  -.093.  As  indicated  in  equation  (28)  and  Table  5,  we 
choose  CONVER  =  .000001. 

g)  Line  10 

As  indicated  in  Figure  9,  X0  is  1.156  ft.  Furthermore , 

Ax  =  .25  inch  (.020833  ft),  thus  DX  =  .020833.  The  throat  height 
of  the  inlet  is  .798  inch  (.0665  ft),  thus  XH1  =  .0665.  The  cowl 
angle  6^  is  7.97  degrees,  thus  DELTAC  =  7.97. 

7 .  S a mpLe  Calculation:  Downstream  Inlet  Region  for  MCAIR  Case  35 

In  Figure  l  ib,  the  coordinate  system  for  the  downstream  inlet  region  of 
MCAIk  Case  35  is  shown  (only  odd-numbered  q-lines  are  shown).  The  upstream 
boundary  of  this  region  corresponds  to  the  i  =  31  station  of  Figure  11a. 

The  input  data  is  indicated  in  Table  7,  and  is  explained  below: 
a)  Line  1 

A  mesh  of  40  points  by  48  points  (M  =  40,  N  =  48)  is  chosen. 

The  maximum  number  of  iterations  permitted  is  MAXITR  =  1000.  Since 
the  downstream  region  overlaps  the  upstream  region  (with  i  =  1  in 
the  downstream  region  chosen  to  be  coincident  with  i  =  31  of  the  up¬ 
stream  region,  i.e.,  IOVLP  =  31),  the  upstream  coordinate  transforma¬ 
tion  data  must  be  read  from  file  OVRLAT  and  thus  1RDZET  =  1.  This 
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TABLE  7.  INPUT  DATA  FOR  COORDINATE  SYSTEM  OF 


DOWNSTREAM  INLET  REGION  (MCAIR  CASE  35) 


Line  1 

M 

N 

MAXITR  IRDZET 

10VLP 

(Format 

615)  40 

48 

1000 

1 

31 

Line  2 

J  H 

J01 

J10 

.Ill 

I  RE  FI 

(Format 

615)  18 

18 

28 

28 

20 

Line  3 

JREF 

IDSCNT 

(Format 

215) 

28 

3 

Line  4 

ID(1) 

ID(  2) 

ID(  3) 

(Format 

1615) 

36 

15 

33 

Line  5 

JD(1) 

JD(  2) 

JD(  3) 

(Format 

1615) 

1 

48 

48 

Line  6 

DY00 

DY01 

DY10 

( Format 

4F10.4 

.0002476 

.00025 

.0002532 

Line  7 

] 

.  .1(1) 

CC2( 1) 

CC1 (M) 

(Format 

4F10.4) 

3.101 

3.984 

2 .83529 

Line  8 

SCALE1 

SCALE 2 

SCALE  5 

(Format 

4F10.4) 

.374 

.374 

1  .5 

Line  9 

S0(1) 

Sl(l) 

CONFER 

( Forma  t 

3F10.4) 

1.78125 

.5  705 

.00O0O1 

Line  10 

X0 

DX 

XU  i 

( Fo  rma  t 

4F10.4) 

1.78125 

.020833 

.066 50 
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I  RE  T 
0 

I  REF 
20 


job  is  the  initial  attempt  to  solve  equations  (22)  in  the  downstream 
region,  and  thus  IRETRY  =  0  (the  equations  (24)  will  converge  in 
this  case  in  less  than  1000  iteration) . 

b)  Line  2 

The  values  of J00  and  J10  (defining  and  at  i  =  1)  are 

chosen  identiacl  to  the  values  used  in  the  upstream  region  (J00  =  18, 

J10  =  28).  A  preliminary  analysis  (not  contained  in  this  report), 

similar  to  that  described  in  Section  III. A. 3,  together  with  program 

BNDRY  indicates  that  J10  =  18  and  Jll  =  28.  Furthermore,  this  analysis 

indicates  DY01  =  DY11  =  .00025  ft,  which  implies  that  y'  is  essentiallv 

m 

unchanged  on  the  ramp  and  cowl  (see  line  6).  Thus,  the  values  of 
IREF1  and  IREF2  are  irrelevant.  They  are  arbitrarily  chosen  to  be 
equal  to  20. 

c)  Line  3 

The  value  JREF  =  28  is  again  found  to  be  satisfactory.  There 
are  three  discontinuities  (see  Figure  lib),  and  thus  IDSCNT  =  3. 

d)  Lines  4,  5 

The  discontinuities  are  located  immediately  downstream  of  the 
following  points:  (i,j)  =  (36,1),  (i,j)  =  (15,48) ,  and  (i,j)  =  (33,48) . 
Thus,  ID(1)  =  36,  JD(1)  =  1,  ID(2)  =  15,  JD(2)  =  48,  ID(3)  =  33  and 
JD(  3)  =  48. 

e)  Line  6 

The  value  of  DY00  is  equal  to  Y^  9  -  Y^  ^  in  the  upstream 
region.  From  the  output  from  the  upstream  region  (see  Figure  12b), 
we  note  that  Y^  2  “  ^33  ^  =  .0002478.  Similarly,  DY01  is  equal  to 
^31  48  ”  ^31  47  w^c*1  from  the  output  is  equal  to  .00025.  A 
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preliminary  analysis  (not  contained  in  this  report)  indicates 

y^  =  2.5  x  10  M  ft.  is  satisfactory  at  i  =  40,  and  thus  DY10  =  DY11  = 

.00025. 

f)  Line  7 

The  values  of  CC1(1)  and  CC2(1)  correspond  to  the  values  of 
Cl  and  C2  at  i  =  31  in  the  upstream  region.  From  the  output  (see 
Figure  12a),  CC1(1)  =  3.101  and  CC2(1)  =  3.084.  An  analysis  of  the 
downstream  boundary  using  BNDRY  (not  contained  in  this  report)  in¬ 
dicates  CC1(M)  =  2.83529  and  CC2(M)  =  2.82212. 

g)  Line  8 

The  values  of  SCALE1  and  SCALE2  are  chosen  as  before.  Since 
y^  is  essentially  uncharged  on  the  ramp  and  cowl,  SCALE3  and  SC.ALE4 
are  irrelevant  and  are  arbitrarily  taken  to  be  equal  to  1.5. 

h)  Line  9 

The  upstream  boundary  (coincident  with  i  =  31  of  the  upstream 
region)  is  located  at  x  =  21.38  inch  =  1.78125  ft.,  thus  S0(1)  = 
1.78125.  At  the  upstream  boundary  on  tie  ramp,  the  mesh  point 
(i,j)  =  (l,48)  is  a  distance  of  6.44  inch  =  .53705  ft.  along  the 
cowl  from  the  cowl  leading  edge  and  thus  Sl(l)  =  .53705.  The  value 
of  CONVER  is  same  as  previously  used. 

i)  Line  10 

The  upstream  boundary  is  located  at  21.38  inch  =  1.78125  ft. 
from  the  ramp  leading  edge,  thus  X0  =  1.78125.  The  remaining  variables 
are  identical  to  their  previous  values.  A  preliminary  analysis 
indicates  Ax  =  .25  inch  =  .020833  ft.  is  satisfactory,  and  thus 
OX  =  .020833. 
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8. 


Sample  Output:  Upstream  Inlet  Region  for  MCAIK  Case  35 


The  output  from  the  program  is  illustrated  in  Figures  12a  through 

12e,  where  sample  of  each  major  section  of  output  for  the  upstream  inlet 

region  of  MCAIR  Case  35  are  given.  In  Figure  12a,  the  values  of  the  input 

data  are  printed,  together  with  the  computed  values  of  CC1(I),  1=1,  ...,  M 

and  CC2(I),  1=1,  . .  . ,  M,  where  the  index  I  increases  from  left  to  right  alr  .ir 

each  row  of  output.  In  Figure  12b,  a  portion  of  the  computed  values  of  the 

coordinates  x..,  y.  .  are  shown.  The  values  of  the  coordinate  pair  (x.  ,,v.  ) 
ij  i,J 

in  feet  are  given  at  each  value  of  n  (i.e.,  at  each  value  of  j)  in  a  sequent'  i. 
manner.  For  each  value  of  n,  the  index  i  increases  from  left  to  right  alone 
each  row  of  output.  The  coordinate  transformation  derivatives  appear  next, 
as  illustrated  in  Figure  12c,  where  the  results  for  £  are  shown.  Again, 

“i ,  j 

for  each  value  of  r|,  the  quantity  £  is  listed,  with  the  index  i  increasing 

Xi,  j 

from  left  to  right  along  each  row  of  output.  Similarly,  the  results  for  , 

rix  and  are  printed.  In  Figure  12d,  the  results  for  XN(I,J)  are  shown 
J  <_  JREF  in  the  same  manner.  These  are  followed  by  X!v(I,J)  for  J  >  JREF  {p.  t 
shown).  Finally,  the  distance  along  the  lower  and  upper  boundaries  (ie.e.. 

S0  and  SI)  are  displayed  as  in  Figure  12e. 

In  calculations  involving  overlapping,  the  message  "R>1  ERROR  0142  on  !.iE 
OVRLAP"  will  appear  ir  the  dayfile.  This  message  is  of  no  significance,  trl 
may  be  ignored. 


SECTION  IV 


CALCULATION  OF  INLET  FLOWFIELD 

A.  Introduction 

The  program  INLET  is  employed  to  compute  the  inlet  flowfield. 

Prior  to  execution  of  this  program,  the  coordinate  transformation  must 
be  calculated  as  discussed  in  Section  III. 

The  following  table  provides  information  on  the  resources  required  to 
execute  program  INLET  on  the  CYBER  175. 

TABLE  8.  RESOURCES  REQUIRED  TOR  EXECUTION 
OF  PROGRAM  INLET 

Resource  Deta_i_Ls 

Computer  time  See  discussion  below 

Input/Output  time  <  60  seconds  (typical) 

Core  Memory  225,000  words  (octal) 

Files  required  INPUT,  OUTPUT,  RSTART, 

STORK,  XY 

(See  Section  IV. F) 

Because  the  total  execution  time  required  to  obtain  a  converged  solution 
within  a  given  mesh  region  is  typically  tens  of  minutes  to  an  hour,  tl  most 
efficient  approach  is  to  submit  a  sequence  of  computer  jobs.  Eacn  job  reads 
in  the  most  recent  results  for  the  flow  variables,  integrates  the  equations 
of  motion  for  a  number  of  time  steps  (see  ITER)  and  then  stores  the  results 
on  file.  litis  procedure  is  repeated  until  a  converged  solution  is  obtained. 
Due  to  the  nature  of  the  code,  the  computer  time  per  time  step  may  vary  con¬ 
siderably  for  a  given  mesh  during  the  course  of  a  calculation.  Typical  values 
are  10  to  20  seconds  per  time  step  on  the  CYBER  175  for  a  high  speed  inlet 
with  a  40  by  48  grid.  The  code  incorporates  a  timer,  which  automa t i ca i 1 v 
determines  at  each  time  step  whether  the  remaining  computer  time  is  sufiicient 
for  continued  execution.  If  sufficient  Lime  does  not  remain,  the  integration 


is  terminated  and  the  results  written  on  file  for  continued  execution  by 
the  next  job.  Thus,  no  job  can  run  out  of  computer  time. 

B.  Description  of  Input  Variables 

The  input  to  program  INLET  consists  of  a  variable  number  of  cards. 
The  definition  and  format  of  the  input  data  is  indicated  in  the  required 
order  in  Table  9.  Additional  information  is  provided  in  Section  IV. D.  Useful 
guidelines  for  the  selection  of  certain  input  parameters  are  provided  in  Table 


TABLE  9.  INPUT  DATA  FOR  PROGRAM  INLET' 


Line  1:  UINF, 

PINF,  XMINF,  CR,  ALPHA,  BETA, 

TMAX 

(Format  7F10.8) 

Fortran 

Definition 

Range  or  Value 

UNIF 

Uco 

Uni ts : 

f  t/sec 

PINF 

Pco 

Uni ts : 

lb  f /ft2 

XMINF 

M 

CO 

CFL 

Courant-Freidrichs- 
Lewy  number 

0<GFL<1 

ALPHA 

a 

0<u-  5 

BETA 

6 

See  Section 

IV. 

TMAX 

CP  time  requested 
on  job  card 

1'n  i  ts : 

seconds 

Line 

2:  XLMBDA(l) , 

XLMBDA( 2) ,  PDOWN,  SREF( 

1) ,  SRL 

F(  2)  | 

(Format:  5FI0.8) 

Fortran 

Definition 

Range  or  Va I 

u_e 

XLMBDA(l) 

X  on  n  =  0 

Units: 

feet 

XLMBDA( 2) 

A  on  n  =  1 

Ua its: 

1  eet 

PDOWN 

static  pressure 
at  C,  =  1 

Un its: 

11)1 /l't  r 

SREF(l) 

s  for  n  =  0 
boundarv  Javer 
(see  (17)) 

-n  i  t  s  : 

1  i’t’t 

SREF(  2) 

s  for  U  =  1 
boundary  layer 
(see  (17)) 

Units: 

feet 

1 

1 

1 


TABLE  9 


CONT '  1) 


Line  3:  IL,  JL, 

ITER,  ILE ( 1) ,  ILE (2) , 

. . — \ 

JSL0,  JSLl  j 

(Format  715) 

Fortran 

Def i n i t  Lon 

Range  or  Value 

IL 

Number  of  points 
in  ((-direction 

<40 

JL 

Number  of  points 
in  redirect  ion 

<48 

ITER 

Number  of  time 
s  teps 

r0 

ILE(l) 

Value  of  i  at 
leading  edge  on 
n  =  0 

1^ILE( 1) <IL 

ILE  (  2) 

Same  as  above  for 

n  =  1 

Same  as  above 

JSL0 

Number  of  points 
in  CSL  on  n  =  0 

3<JSL(0<2O 

JSL1 

Same  as  above  for 

n  =  l 

3;- .  I SL 1;  20 

Line  4:  JREF0,  JREF1 ,  I DAMP 1,  I DAMP 2,  I START,  TPRINT 

(Format  615) 


Fortran 

Def ini t ion 

Range  or  Value 

JREF0 

See 

Section 

IV.  D 

JREF1 

See 

Sec  t i on 

IV.  D 

I DAMP 1 

See 

Sect i on 

IV.D 

0  o  r  i 

I DAMP 2 

See 

Sect i on 

IV.  I) 

0  or  1 

I START 

See 

Section 

IV.D 

0,  J  or  2 

I PR I NT 

Flow  data  printed 
every  IPR1NT  steps 

1 
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TABLE  9.  CONT'D 


Line  5:  ITNS(l),  ITNS(2),  ITNE(l),  TTNE(2),  IOVLP,  IP DOWN 

(Format  61  ,) 


Fortran 

Def ini ti on 

Rang 

e  or  Value 

ITNS(l) 

Value  of  i  at  start 
of  transition  on 

n  =  o 

l  < !  'j : 

SS(i);.IL 

ITNS(2) 

Same  as  above  for 

n  =  1 

Sana1 

as  above 

ITNE(l) 

Value  of  i  at  end 
of  transition  on 

n  =  o 

Same 

as  abi  ve 

I1’NE(  2) 

Same  as  above  for 

n  =  1 

Sap, 

:!:■>  aho  «; 

IOVLP 

See  Section  IV. D 

IP DOWN 

See  Section  IV. D 

0  or 

I 

(Format  7F10.8) 

Fortran  jlcfinition  K.un;.  >u  V  a  hie 

XMDOT (1,1)  ,  m  on  n  =  0  I'nu-  o  sluy-.s/lt 

I  =  1 , . .  . , IL 

XMDOT  (1,2),  m  on  v\  -  1  Same  as  .ib  -ve 

1  =  1,  .  .  .  ,  IL 


TABLE  10.  GUIDELINES  LOR  INPUT  DATA  FOR  PROGRAM  INLET 


Parameter (s ) 

UINF ,  PINF,  XMINF 

CFL 

ALPHA 

BETA 

XLMBDA(l),  XLMBDA(2) 

PDOWN 

SREF(l),  SREF (2) 

ILE(K) 

JSL0,  JSL1 
JREF0,  JREF1 
ITNS(K),  ITNE(K) 
IOVLP 
IPDOWN 
XMDOT(I ,K) 


Gui  delaine 

These  quantities  are  the  f roes t ream  values 
ahead  of  the  inlet  entrance. 


Typically,  a  value  of  0.9  is  satisfactory. 
Occurrence  of  numerical  instability  (manifested 
by  arithmetic  overflow  error  messages)  can  often 
be  cured  by  restarting  the  calculation  and  using 
a  smaller  value  of  CFL,  although  code  efficiency 
is  thereby  decreased. 

The  pressure-damping  term  (55),  (56)  prevents 
numerical  instability  in  the  presence  of  shock 
waves  and  eliminates  unphysical  "wiggles”  in 
the  pressure  therein.  A  value  of  a  =  5.0  is 
typical.  Too  high  a  value  of  a  can  cause 
numerical  instability. 


Must  be  used  only  to  control  transients  if 
necessary.  Over  at  least  several  t  (see  (57)), 
3  must  equal  1. 

When  experimental  values  for  ,S0  in  (17)  are  not 
available,  the  calculation  can  be.  run  awhile 
with  no  relaxation  until  a  good  estimate  for 
<50  is  obtained.  If  relaxation  is  not  used, 
set  both  parameters  equal  to  1. 

See  Note  5  in  Section  IV. 1).  If  fPlXtWN  -  0, 
set  PDOWN  =  0. 


For  no  relaxation,  set  SREF(l)  and  SREF(2) 
equal  to  anything  greater  than  S(IL,1)  and  S(TL,2), 
respectively  (e.g  ,  999).  For  relaxation,  set 
SREF  equal  to  arc  distance  at  which  relaxation 

begins . 

See  Note  1  in  Section  1V.D. 

See  Note  4  in  Section  IV.;). 


See  Note  3  in  Section  IV.!'). 


See  Note  2  in  Section  iY.D. 


If  not  overlapping,  set  IOVLP  =  0. 
See  Note  5  in  Section  IV.  I). 

See  Note  7  in  Section  IV. I). 
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C.  Flow  Charts 


1 .  Overall  Program 

The  pertinent  features  of  the  overall  program  INLET  are  indicated 
below.  The  names  of  the  subroutines  called  are  indicated  in  capital 
letters  in  parentheses.  CSL  refers  to  computational  sublayer. 


Subroutine 

INIT 


BEGIN 


Update  CSL  on  n  =  0  and 
n  =  1  (VSL) 

Update  boundary 

conditions  (BC) 

ITER  =  0?  u 


Yes 


1 


No 


BEGIN  ITERATION  LOOP 


DO  1200  NADV  =  1,  ITER 


2.  L  Operator  (subroutine  LZETA) 

For  each  region,  the  L  operator  is  applied  to  a  domain  defined  by 
2_f_i_f_IL-l  =  12,  JS  <  j  <  JE .  The  pertinent  features  are  indicated 
below • 


BEGIN 


PREDICTOR  STEP 


Compute  F  .  (F) 
-L  9  J 

_* 

Compute  F.  .  (F) 
1  >  3 


Form  U. 


i.J 


DO  for 
i  =  2,  12 


DO  for 
j  =  JS,  JE 


Compute  u,  v,  e^,  p 
Apply  boundary  conditions  (BC) 
Update  e  (EPSLN) 


CORRECTOR  STEP 


Compute  F  .  (F) 
^  *  J 


Compute  F  (F) 

»  .J 


Form  U 


i .  i 


DO  for 
i  =  2,  12 


DO  for 
.1  =  JS,  JE 


Compute  u,  v,  e^,  p 
Apply  boundary  conditions  (BC) 
Update  e  (EDDY) 


RETURN 
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3. 


^Operator  for  Regions  1  and  2  (subroutine  LETAL) 

The  £■  operator  is  applied  to  a  domain  given  by  2  <_  i  ^  IL  -  1  =  12, 
JS  j  <_  JE.  The  pertinent  features  of  the  operator  for  Regions  1  and 
2  are  indicated  below.  The  1.  operator  for  the  other  regions  is  very 
similar  and  consequently  is  not  detailed  separately. 


BEGIN 


DO  for 

=  2,  12 


DO  for 
j  =  JS,  JE 


PREDICTOR  STEP 

-- i-n 

.Is  IRGN  =  1?J - No 


Yes 


Compute  G.  ( G ) 

1,1 


Compute  j  (G) 


** 

Fo  rm  U . 

i,  j 


CORRECTOR  STEP 


1 


,  IRGN  =  1?  r 

_r_ 


Yes 


Compute  G.  (G) 
1 ,  z 


** 

Compute  G^  Jc, 
(GCON) 

DO  for 
=  JS,  JE 


** 


Compute  G.  (G) 

1 » J+J- 


** 
1.  . 
i,3 


Form  U. 


RETURN 


D.  Source  Code  Notation 


The  definition  of  the  pertinent  source  code  variables  is  given  below. 
For  any  array,  the  index  I  is  equivalent  to  the  subscript  i,  and  the  index 
J  is  equivalent  to  the  subscript  j.  The  index  K  assumes  values  of  1  or  2, 
with  K  =  1  referring  to  the  n  =  0  boundary  layer  and  K  =  2  referring  to 
the  n  =  1  boundary  layer.  The  abbreviation  CSL  refers  to  the  computational 
sublayer,  and  OM  refers  to  ordinary  mesh. 


Fortran 

Definition 

Range  or  Value 

Reference 

ALPHA 

a  (pressure  damping) 

0  <  a  <  5 

(55),  (56) 

BETA 

6 

Q  <  0  (transients 

** 

1 1 .  B  .  5 

only) 

6  =  1  (converged 

sol ution) 

CAPGAM(l) 

,  CAPGAM0 

A  for  n  =  0 

Units:  feet 

(16) 

CAPCAM(2) 

,  CAPGAM1 

A  for  n  =  1 

Units:  feet 

(16) 

CF 

,1  „2 
=  t  !-=  p  U 

I  W  z  oo  oo 

CFL 

Courant-Friedrichs- 

0  <  CFL  ••  1 

(44),  (45) 

Lewy  number 

TT .B  .1 

CV 

c^  (air) 

4290  ft2/see2-°R 

I  .C  .2 

DETA 

An 

(21) 

DTI,  DT2, 

DT3, 

Quantities  employed  in 

DT4,  DT5 

determining  proper  time 

step  in  each  region 

DT 

At 

Units:  seconds 

(37) 

DUDN 

3a/ 3n 

(ID 

*** 

DY  PLUS 

y  '  u .  /  v 

(61) 

DY SL  PLUS 

*** 

Ay  u.N/v 

SL  *  w 

(14),  TIT . 

T.e.,  equations  (55)  and  (56).  T.e.,  Section  1T.R.5. 

*** 

Appear  on  program  output. 


Fortran 

Definition 

Range  or 

Value 

Reference 

DZETA  A? 

(21) 

EFFIC 

Measure  of  efficiency 

(51) 

of  code 

2  , 

2 

EI(I,J) 

e  at  (x  ,y  ) 

y  J  x  >  J 

Units : 

ft  /sec 

I.C.2 

EINF 

e 

00 

Units: 

,  2, 
ft  /sec 

2 

I.C.2 

EIINF 

eia> 

Units : 

,  2, 
ft  /sec 

2 

1  .C .  2 

EPSLN(I, J) 

z  ,  .  in  OM 
i  >  J 

Units : 

lbf  sec 

/ft2 

(15),  (17) 

EPSLNI 

c  ^ 

Units: 

lbf  sec 

/ft2 

(11) 

e<li 

EPSLNO 

E  e<^o 

Units : 

lbf  sec 

/ft2 

(12) 

EPSLSL(I,J,K) 

e .  .  in  CSL 
i>  j 

Units: 

lbf  sec 

/ft2 

(15),  (17) 

EREF(J) 

e  (s  ,n)  in  OM 

Units : 

lbf  sec/ft 

(17) 

eq  o 

ESLREF (J,K) 

e  (s  ,n)  in  CSL 
eq  o 

Uni ts : 

lbf  sec 

/ft2 

(17) 

ETAX(I,J) 

nxiJ 

Units : 

ft"1 

I.C.2 

ETAY(I,J) 

nyi,j 

Units : 

ft'1 

I.C.2 

FAC (I ,K) 

A  r  [  p  /  26p 

Units : 

ft"1 

1  w1  w  w 

GAM(I,K) 

r 

(16) 

GAMMA 

Y 

1,4  (air) 

(10) 

j 

0:  no  pressure 

IDAMP1  < 

damping  for  d 

1  _  _  _  _  ^ 

(55),  (56) 

1:  pressure 

damping  for 

IDAMP2 

Same  as  above  for  L 

n 

IL 

Number  of  points  in 

£40 

4-direction 

12 

IL  -  1 

ILE(l) ,  ILE0 

Value  of  i  at 

1  £  II.E0 

£  IL 

* 

Note  1 

leading  edge  on  n  = 0 

ILE(2),  ILE1 

Same  as  above  for 

1  £  ILE1 

£  IL 

Note  1 

i 

1  =  1 

IOVLP  Value  of  i  in  upstream 

1  <  IOVLP  <  IL 

-  5 

II.  B.  3 

region  that  defines 
restart  station 


Notes  are  listed  at  end  of  this  section. 
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Fortran 


Definition 


IPDOWN 


IPRINT 


Range  or  Value  Referen 


0: 

no  downstream 

Note  S 

pressure  applied 

1: 

pressure  PDOWN 

applied  at  C  =  1 

Flow  data  printed 
every  IPRINT  steps 


IRGN 


I START 


Region  in  split  mesh  1  to  a 

r 

0:  Initialize  vari¬ 

ables  to  uniform 
flow  (file  RSTART 
not  used) ;  i .e . , 

Uoo,  Poo,  etc. 

1:  Continue  calcula¬ 

tion  using  data 
<  from  file  RSTART 

2:  Read  upstream  region 

from  file  RSTART. 

Set  i  =  IOVLP  as 
restart  station  of 
downstream  region 
and  commence 
calculation 


II  .B  .2 


ITER 

ITNF(l),  ITNE0 

TTNE(2),  TTNF.l 
ITNS(l),  ITNS0 

ITNS (2) ,  ITNS1 

JI1 


Number  of  time  steps 
(if  ITER=0,  code  simply 
prints  results  with  no 
integration ) 

Value  of  l  at  end  of 
.  .  .  * 
transition  region 

on  n  =  0 

1  <  ITNS0 

<  1TNE0  <  IT. 

Note 

Same  as  above  for 

n  =  l 

1  <_  1TNS1 

<  ITNE.1  <  II. 

Note  2 

Value  of  i  at  begin¬ 
ning  of  transition 
region*  on  n  =  0 

See  1TN1.0 

ah  live 

Note  2 

Same  as  above  for 
n  =  1 

See  ITNE1 

above 

Note  2 

Upper  limit  for  j 
in  Region  1 

T1  .B  .2 

* 

Region  of  transition  from  laminar  to  turbulent  flow. 
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Fortran 

Definition 

Range  or  Value 

Referenc 

JI2 

Upper  limit  for  j 
in  Region  2 

II. B. 2 

JI3 

Upper  limit  for  j 
in  Region  3 

II  .B .  2 

JI4 

Upper  limit  for  j 
in  Region  4 

II. B. 2 

JL 

Number  of  points  in 
n-direction 

9  £  JL  <  48 

J2 

JL  -  1 

JMATCH (1,1) 

Value  of  j  in  OM  at 
which  e  transfers  from 
inner  to  outer 
expression  for  n  =  0 
boundary  layer 

2  to  JREF0 

(15) 

JMATCH(I , 2) 

Same  as  above  for 
n  =  1  boundary  layer 

JREF1  to  J2 

(15) 

JREF0 

For  all  j  £  JREF0, 
the  eddy  viscosity 
refers  to  the  n  =  0 
boundary  layer 

Note  3 

JREF1 

For  all  j  _>  JREF1 , 
the  eddy  viscosity 
refers  to  the  n  =  1 
boundary  layer 

No  te  3 

JSL0 

Number  of  points  in 

CSL  on  n  =  0 

3  £  JSL0  £  20 

No  te  4 

JSL1 

Same  as  above  for 
n  =  1 

3  <  JSL1  <  20 

Note  4 

M(l) ,  Ml 

ml 

>1 

(47) 

M(2) ,  M2 

m2 

>1 

(47) 

M(  3) 

m3 

1 

(47) 

M(4),  M4 

m4 

>1 

(47) 

M(  5) ,  M5 

m5 

>1 

(47) 

NADV 

Time  step  number 

85 


Fortran 

Definition 

Range  or 

'  Value 

Reference 

NI 

Cumulative  number  of 
time  steps 

P(I,J) 

Static  pressure  p  . 

i  >  3 

Units : 

lbf/ft2 

I  .C  .2 

PDOWN 

Static  pressure  at  c,  =  1 

Uni ts : 

lbf/ft2 

Note  5 

PINF 

Poo 

Units : 

lbf/ft2 

PR 

Pr 

0.72  (air) 

(10) 

PRTURB 

Prt 

0.90 

(10) 

QYWL ( I , K) 

Heat  transfer  at 

0  (adiabatic  wall) 

wall  in  CSL 

QY2(I,K) 

Heat  transfer  normal 
to  wall  at  matching 
point  in  CSL 

Units: 

Ibf/ f t-sec 

II  .8.7 

REM(l) 

Max  Re •  in  CSL  on  n  =  0 
m 

Must  be 

less  than 

(65),  (66 

approximately  0.25 

REM(  2) 

Same  as  above  for  CSL 
on  n  =  1 

Same 

(6  5)  ,  l6r 

RHO(I, J,L)* 

Pi,j 

Units: 

sl  ugs/  ft 2 

(3)  to  ... 

RHOE(I,j,l)* 

P6i,  j 

Units : 

lbf/ft2 

(3)  tv'  (.•' 

RHOINF 

Pco 

Units : 

S.1  Ugs  /  ft2 

RHOU(I, J,L)* 

pu.  . 

1>J 

Un its: 

2 

slugs/ f t^-sec 

(  3)  to  l  C 

RHOV ( I , J,L) * 

PV. 

1.  j 

Units: 

O 

s] ugs/f t  -sec 

(3)  to  (o 

RK1 

K 

0.40 

(11) 

RK2. 

k2 

0.0168 

02) 

S(I,1) 

Distance  along  n  =  0 

Units : 

fee  t 

See  S0(V: 

from  leading  edge 

in  11  1  . ( '  . 

S(I,2) 

Distance  along  t,  =  1 

Un i ts: 

f  ee  t 

See  SIU' 

from  leading  edge 

in  Ti  l  .(' 

* 

L  =  1 

represents  value  after  corrector 

step,  and 

I.  =  7  represents 

value  after 

predictor  step.  For  any  analysis 

of  conve 

rged  solution. 

always  use  L 

=  1  level . 

Fortran 

Definition 

Range  or  Value 

Referent 

SREF(l) 

s  for  n  =  0 

0 

boundary  layer 

Units : 

feet 

(17) 

SREF (2) 

Same  as  above  for 
n  =  1  boundary  layer 

Units: 

feet 

(17) 

TAUWL (1,1) 

Wall  shear  stress  on 
n  =  0 .  Positive  values 
imply  force  is  active  in 
positive  C-direction, 
and  vice  versa 

ITNS0  < 
Units : 

i  <  IL 
lbf/  ft  2 

TAUWL (1,2) 

Same  as  above  for  n  =  1 

ITNS1  < 
Units : 

i  <  IL 
lbf/ft2 

TAU2(I, 1) 

Shear  stress  parallel 
to  wall  in  CSL  on  n  =0 
at  matching  point 

Same  as 

TAUWL (1,1) 

TAU2(I,2) 

Same  as  above  for  n  = 1 

Same  as 

TAUWL (I, 2) 

TEMP 

Static  temperature 

Degrees 

Rankine 

TIME 

Cumulative  physical 
time  of  integration 

Units : 

seconds 

TMAX 

CP  time  requested  on 
job  card 

Units : 

s  or ends 

TSL(I,J,K) 

Static  temperature 
in  CSL 

Units: 

Degrees  Rankine 

U(I,J) 

u.  .  in  OM 

i.j 

k 

Uni ts: 

ft/ sec 

(A) 

UDELTA,  UREFDL 

U  6. 

rer  1 

Uni ts : 

ft “/sec 

(12) 

UINF 

Uo,  =  magnitude  of 

freestream  velocity 

Units : 

f  t  /  s  ee 

UREF,  UREFER 

Lref 

Units  : 

f  t/ser 

(12) 

USL ( I , J , 1 ) 

u  in  CSL  on  n  =  0 

1TNS0  ■ 
Un i ts : 

i  IL 

ft/ sec 

(18) 

lTSL(  I ,  J,  2) 

u  in  CSL  on  ri  =  1 

ITNS1  < 
Uni ts : 

1  <  Tl. 
ft /sec 

V(I,.J) 

v.  .in  OM 

1.1 

Units : 

ft /sec 

(A) 

X(I,J) 

x.  .in  OM 
i,.l 

Uni ts : 

feet 

r 


Fortran 

Definition 

Range 

or  Value 

Reference 

XJINV(I, J) 

1/J.  . 

1,3 

Units: 

feet 

(7) 

XLMBDA(l) ,  XLMBDA0 

A  on  n  =  0 

Units : 

feet 

(17) 

XLMBDA(2),  XLMBDA1 

A  on  q  =1 

Units : 

feet 

(17) 

XMDOT(I, 1) 

m^  on  q  =0 

Units : 

slugs/ f t^-sec 

Note  7 

XMDOT (I , 2) 

a,  on  n  =  1 

Units : 

slugs/ f  t^-sec 

Note  7 

XMINF 

Freestream  Mach  No 

.  M 

co 

XMUSL (I , J, 1) 

y  in  CSL  on  q  =  0 

ITNS0 

<  I  IL 

(8) 

Units : 
(air) 

lbf  sec/ ft 

XMUSL (I , J ,  2) 

Same  as  above  for 

n  =  1 

ITNS1 

<  I  <  IL 

(8) 

Units: 

(air) 

lbf  sec/f t 

XN ( I , J ) 

See  XN ( I , J )  in  III.C.4 

Y(I,J) 

y.  .  in  OM 
i,j 

Units: 

feet 

ZETAX(I, J) 

Units: 

ft"1 

ft-1 

I.C.2 

ZETAY(I, J) 

?yi,3 

Units : 

I .  C .  2 

Note  1:  On  q  =  0, 

boundary  conditi <np 

(53) 

for  a  nc 

>-sl ip  wal 1  are 

applied 

for  i  >_  ILE(l)  .  If  the  entire  surface  n  =  0  is  a  no-slip  surface, 
ILE(l)  =1.  If  the  entire  surface  is  within  the  freestream,  set 
ILE(l)  =  IL .  Similar  results  hold  for  n  =  1. 


Note  2:  If  the  entire  boundary  layer  on  q  =  0  is  fully  turbulent,  set 
ITNS(l)  =  ITNE(l)  =1.  If  the  entire  boundary  layer  on  ri  =  0 
is  fully  laminar  or  if  there  is  no  boundary  layer,  set  ITNS(l)  =  IL 
and  ITNE(l)  =  IL .  Note  that  in  mesh  overlapping  a  restart  station 
cannot  fall  within  the  region  of  transition  from  laminar  to  turbulent 
flow  in  any  boundary  layer.  Similar  results  hold  for  q  =  1 . 

Note  3:  In  choosing  JREF0  and  JREF1,  the  following  criteria  must  be  observed: 

a)  1  <  JREF0  £  JREF  <  JREF1  <  JL ,  where  JREF  is  defined  in  Section 

III.C.2 

b)  The  contour  n  =  HjreF0  roust  lie  outside  of  the  boundary  layer  on 
ri  =  0  at  all  stations.  The  contour  need  not  be  close  to  the 
actual  boundary  layer  thickness  at  all .  For  boundary  layers  on 
the  ramp  and  cowl  of  roughly  equal  thickness  at  any  station,  .TREE 
is  chosen  such  that  HJREF  roughly  bisects  the  inlet  height  at  all 
stations.  This  can  be  estimated  using  the  mesh  distribution  on 

C  =  0  and  C  =  1 . 

c)  The  contour  q  =  PjreFI  roust  lie  outside  the  boundary  layer  on 

n  =  1 . 


Note  4: 


Note  5 


Note  6 


Note  7 


1 


The  quantities  JSL0  and  JSL1  cannot  be  changed  during  the  course 
of  any  calculation,  nor  during  overlapping.  Since  the  CSL 
calculation  is  very  fast,  it  is  recommended  to  use  JSL0  =  JSL1  =  20 
in  all  cases. 

The  effect  of  a  terminal  shock  can  be  achieved  by  specifying  a 
desired  downstream  pressure  PDOWN  (lbf/ft^)  and  setting  IPDOWN  =  !. 
For  sufficiently  high  downstream  pressures,  a  terminal  shock  will 
form  in  the  inlet  throat.  The  necessary  value  of  PDOWN  can  be 
approximated  by  a  simple  normal  shock  analysis  using  data  at  r_  -  1 
(when  solution  is  converged  and  IPDOWN  =  0),  but  the  location  and 
stability  of  the  terminal  shock  will,  in  general,  be  somewhat 
sensitive  to  the  value  of  PDOWN,  as  in  the  experimental  case. 

:  For  all  CSL  variables  (e.g.,  USL (I ,  J,K.) )  ,  the  index  I  refers  to 

the  location  in  the  4-direction  in  the  same  manner  as  all  other 
variables  (e.g.,  U(I,J)).  The  J  index  is  defined  such  that  J  =  1 
is  always  the  surface  (i.e.,  n  =  0  or  n  =  1)  and  J  =  JSL0  or 
JSL1  is  the  matching  point. 

:  The  array  XMD0T(I,K)  (i.e.,  mi)  is  the  bleed  mass  flux  (in  units 

of  slugs/f t^-sec)  at  the  station  i.  It  is  defined  such  that 
negative  values  always  imply  bleed  (for  either  n  =  0  or  n  =  1 
surface),  and  positive  values  always  imply  blowing.  The  bleed 
mass  flux  is  assumed  uniform  within  each  bleed  zone  and  its  value 
is  given  by 

(Bleed  flow  rate  in  slugs/sec) 

m  =  - - — - - - - — - 2 — - - 

Area  in  ft 

with  proper  consideration  for  the  sign  as  indicated  above.  Note 
that  the  Area  above  is  the  total  surface  area  of  the  bleed  zone, 
not  the  surface  area  of  the  bleed  holes. 

In  executing  the  program,  the  values  of  m^  may  be  changed  as  desired 
to  investigate  effects  of  different  bleed  distributions  on  the  flow 
structure.  Changes  should  be  made  gradually;  e.g.,  by  a  sequence  of 
computer  runs  in  which  is  changed  by  10%  or  so  and  then  the  flow 
is  integrated  for  10  steps,  until  the  desired  mj  is  achieved. 

Note  that  all  values  of  XMD0T(I,K),  I  =  1,IL,  must  be  specified, 
including  those  positions  where  it  is  zero. 
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E.  File  Structure 


The  program  employs,  in  general,  three  files  in  addition  to  the 
conventional  files  INPUT  and  OUTPUT.  The  descriptions  are  indicated  below. 


File 

Description 

Data  Structure 

RSTART 

Input  file  with  latest  results  for 
flow  variables,  or  converged  results 
of  upstream  region  in  the  case 
of  overlapping 

RHO,  RHOU,  RHOV,  RHOE,  U, 
V,  El,  P,  EPSLN,  USL,  TSI. 
EPSLSL ,  XML'S L,  TAUWL,  TAP 
QYWL,  QY2,  FAC,  TIME,  Nl, 
EREF,  ESLREF 

STORE 

Output  file  for  flow  variables 

Same  as  above 

XY 

Input  file  for  coordinate 
transformation  data 

X,  Y,  ZETAX,  ZETAY,  ETAX, 
ETAY,  XN,  S 

The  data  are  written  unformatted  in  the  sequence  indicated.  The 
following  indicates  the  size  of  each  array,  where  for  example  RHO  is 
dimensioned  as  RHO(40,48,2)  . 


Array 

RHO,  RHOU,  RHOV,  RHOE 

Size 

40  x  4  8 

x  2 

U,  V,  El,  P,  EPSLN 

40 

X 

00 

USL,  TSL,  EPSLSL,  XMUSL 

40 

X  20 

x  2 

TAUWL,  TAU2,  QYWL,  QY2,  FAC 

40 

x  2 

TIME,  NI 

1 

EREF 

48 

ESLREF 

20 

x  2 

For  example,  the  sequence  in  which  data  are 

written  on 

f  i  1  e 

RSTART  is 

(consult  any  Fortran  IV  manual) 


KHO(l,l,l), 

RH0(2,1,1) , 

...,  RHO(40,1,1) 

RH0(1 , 2, 1) , 

RH0(2, 2,1) , 

...,  RHO(40 , 2, 1) 

RH0(1,1, 2)  , 

RH0(2,1, 2) , 

i  •  •  • 

....  RHO(40,1,2) 

so  forth. 

The  utilization 

of  the  files 

is  indicated  as  follows: 

File 

Type 

Utilized 

RSTART 

Input 

For  ISTART  =  1  or  2  only 

STORE 

Output 

Required  if  ITER  >  0 

XY 

Input 

Required  for  every  code 
execution 

F.  Sample  Calculation:  Upstream  Inlet  Region  for  MCAIR  Case  35 

In  Figure  9,  the  general  physical  features  of  the  upstream  inlet  region 
of  MCAIR  Case  35  are  indicated.  A  coordinate  system  for  this  region  was  com¬ 
puted  in  the  manner  discussed  in  Section  III.C.6,  except  that  the  mesh  spacing 
Ax  was  taken  equal  to  0.125  inch  and  the  upstream  boundary  at  x  =  13.93  inch. 
The  downstream  boundary  is  thus  x  =  18.81  inch.  Within  this  region,  there- 


fore,  the  boundary  layer  bleed 

schedule 

is 

as 

follows  (See  Figure  9)  : 

Surface 

Extent 

Range 

in 

JL 

m( slugs/ f  t2-sec) 

Ramp 

13 .9  3<x<16 .  5 

i  =  1 

to 

21 

0. 

Ramp 

16 . 5<x<18.8 

i  =  22 

:  to  40 

-3.85  x  10  3 

Cowl 

13.93<x<18.8 

i  =  1 

to 

40 

0. 

The  flowfield  calculation  was  performed  by  submitting  a  sequence  of  several 
jobs  requesting  900  seconds  each  and  having  the  input  data  indicated  in  Table 
11.  Each  input  line  is  described  below, 
a)  Line  1 

The  freestream  velocity  is  2230  ft/sec  (U1NF  =  2230), 
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TABLE  11.  INPUT  DATA  FOR  FLOWFIELD  CALCULATION  OF 


UPSTREAM  INLET  REGION  (MCAIR  CASE  35  WITH 
Ax  =  0.125  INCH) 


Line  1 

UINF  PINF 

XMINF 

CFL 

ALPHA  BETA  TMAX 

(Format 

7F10.8) 

2230.  91 

.84 

3.51 

0.9 

5.0 

1.0  900. 

Line  2 

XLMBDA(l) 

XLMBDA(2) 

PDOWN 

SREF(l) 

SREF(2) 

(Format 

5F10.8) 

1.0 

1 

.0 

0. 

999  . 

999  . 

Line  3 

IL  JL 

ITER 

ILE(l) 

ILE(2) 

JSL0  JSL 1 

(Format 

715) 

40  48 

40 

1 

10 

20  20 

Line  4 

1 

JREF0  JREF1 

I DAMP 1 

I  DAMP  2 

ISTART  IPRINT 

(Format 

615) 

28 

29 

1 

1 

1. 

40 

Line  5 

I 

ITNS(l) 

ITNS (2) 

ITNE(l) 

LTNE(  2)  I.OV’LP 

IP DOWN 

(Format  615) 

1 

20 

1 

30 

I 

0 

Lines  6 

to  11 

XMDOT (1,1) , 

I  = 

1,  ... 

,  40 

(Format 

7F10.8) 

0 

0. 

0. 

0. 

0. 

0. 

0  . 

0 

0. 

0. 

0. 

0. 

0  . 

0  . 

0 

0. 

0. 

0. 

0. 

0. 

u. 

- 

00385  — - 

— 

- — - 

-  -.00  <85 

- 

00385  — - 

—  -.00385 

Lines  12  to  17 

XMD0T(I,2) , 

I  =  1 

,40 

(Format 

7F10.8) 

0. 

0. 

0. 

0. 

0. 

0. 

0  . 

0. 

0  . 

0. 

0. 

0. 

0. 


-  0. 
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Mach  number  is  3.51  (XMINF  =  3.51),  and  pressure  p^  is  91.84 
lbf/ft2  (PINF  =  91.84).  The  Courant  number  is  0.9  (CFL  =  0.9),  and 
the  damping  coefficients  a  and  8  are  taken  to  be  5.0  and  1.0.  The 
requested  computer  time  per  job  is  900  sec  (TMAX  =  900). 

b)  Line  2 

As  discussed  in  Ref.  18,  the  eddy  viscosity  relaxation  model  was 
not  employed  in  the  presence  of  bleed.  This  is  accomplished  by 
setting  SREF(l)  and  SREF(2)  to  any  number  greater  than  the  largest 
value  of  S(I,1)  and  S(I,2).  In  this  case,  SREF(l)  =  SR£F(2)  =  999. 
is  sufficient.  Since  XLMBDA(l)  and  XLMBDA(2)  are  thus  irrelevant, 
they  are  set  equal  to  one.  A  downstream  pressure  is  not  imposed 
(see  Line  5) ,  thus  PDOWN  is  irrelevant  and  arbitrarily  set  equal  to 
zero. 

c)  Line  3 

A  mesh  of  40  points  by  48  points  is  used  (IL  =  40,  JL  =  48) .  A 
total  of  40  time  steps  are  requested  in  each  job  (ITER  =  40) .  The 
leading  edge  of  the  ramp  lies  upstream  of  the  upstream  boundary, 
thus,  ILE(l)  =  1.  The  leading  edge  of  the  cowl  lies  at  x  =  15  inch . 

For  the  mesh  spacing  used,  this  point  lies  between  i  =  9  and  i  =  10 
on  the  upper  surface,  and  thus  ILE(2)  =  10.  A  total  of  twenty  points 
were  used  in  each  CSL  (JSL0  =  JSL1  =  20)  . 

d)  Line  4 

In  the  coordinate  system  generated  for  this  case,  the  value 
JREF  =  28  was  used.  Thus  we  choose  JREF^  =  28  and  JREF1  =  29  (see 
Note  3  of  Section  IV. D.) .  The  parameters  IDAMP1  =  IDAMP2  =  1  indicates 
that  pressure  damping  is  employed  due  to  the  presence  of  shock  waves 
in  the  flow.  The  pai imeter  ISTART  =  1  for  all  runs.  For  the  first  job. 
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file  RSTART  isthe  same  as  file  FLOWDN  from  program  UPSTRM  discussed 
in  Section  IV. H.  For  all  succeeding  jobs,  file  RSTART  is  the  STORE 


file  from  the  previous  job.  We  choose  to  print  the  results  after 
forty  time  steps  (IPFr.NT  =  4C)  . 

e)  Line  5 

Since  the  boundary  layer  on  the  0=0  boundary  is  fully  turbulent 

at  the  upstream  boundary,  ITNS(l)  =  ITNE(l)  =  1  (see  Note  2  of  Section 

14 

IV. D).  Using  the  method  of  Deem  and  Murphy'  ,  it  is  estimated  that 
the  region  of  transition  from  laminar  to  turbulent  flow  on  the  n  =  1 
surface  corresponds  to  the  range  i  =  20  to  i  =30,  and  thus  ITNS(2)  = 
20  and  ITNE(2)  =  30.  Since  ISTART  =  1,  the  value  of  I0VLP  is  irrele¬ 
vant  and  thus  IOVLP  is  arbitrarily  set  to  one.  A  downstream  pressure 
is  not  applied  since  the  terminal  shock  is  assumed  to  be  downstream  of 
the  region,  and  thus  IPDOWN  =  0. 

f)  Lines  6  to  11 

Using  the  bleed  schedule  indicated  above,  the  bleed  mass  flux 
XMD0T(I,1),  1=1,  . . . ,  40  on  n  =  0  is  entered. 

g)  Lines  12  to  17 

Same  as  above  for  T|  =  1 . 

G.  Sample  Output:  Upstream  Inlet  Region  for  MCAIR  Case  11 

The  output  from  the  program  is  illustrated  in  Figures  i 3a  through  1 3f , 
where  samples  of  each  major  section  of  output  for  the  upstream  inlet  region 
of  MCAIR  Case  35  are  presented  sequentially.  The  output  was  generated  by 
one  of  the  later  jobs  in  the  job  sequence.  For  purposes  of  brevity,  ITER 
and  IPRINT  were  changed  to  10. 

In  Figure  13a,  the  values  of  the  input  parameters  are  indicated.  They 
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are  followed  by  a  list  of  variables  printed  at  every  time  step.  These 
variables  are  the  time  step  number  (NADV) ,  the  total  elapsed  physical  time 
(sec)  of  calculation  ror  the  region  (TIME),  the  time  step  (DT) ,  the  values 
of  in  equation  (47)  ,  the  values  of  j  for  the  upper  limit  of  each  of  the 
first  four  mesh-split  regions  (JI1  to  JI4)  ,  the  time  steps  in  each  of  the 
five  mesh-split  regions  (DTI  to  DT5)  and  the  parameter  EFFIC  discussed  in 
equation  (51).  As  indicated  previously,  the  time  step  DT  may  vary  substan¬ 
tially  during  program  execution,  although  the  term  EFFIC  remains  approxi¬ 
mately  constant.  At  the  end  of  the  specified  number  of  iterations,  the 
physical  time  in  seconds  (TIME)  is  specified,  which  is  important  to  the  de¬ 
termination  of  convergence  as  discussed  previously. 

In  Figure  13b,  the  flow  variables  on  the  ordinary  mesh  are  given  at  each 
value  of  i  (that  is,  at  each  station)  denoted  by  "COLUMN  1",  "COLUMN  2",  etc. 
For  brevity,  only  those  results  for  i  =  1  and  i  =  2  are  shown.  At  each  value 
of  i,  the  corresponding  value  of  x  is  indicated  in  feet  ("X  EQUALS  ..."). 

The  values  of  the  cartesian  y  coordinate  (Y)  are  listed  at  each  value  of  j, 
together  with  the  values  of  the  static  temperature  in  °R  (TEMP),  cartesian 
velocity  components  v  and  u  (V,U)  in  ft/sec,  static  pressure  (P)  in  lbf/ft2 
and  turbulent  eddy  viscosity  (EPSLN)  in  lbf-sec/ft2. 

In  Figure  13c,  the  values  of  Urgf  (UREFER)  and  6^  (UDELTA)  are 

listed  at  each  station  for  the  lower  (ETA  =  0)  and  upper  (ETA  =  1)  boundary 
layers.  The  values  of  j  at  which  the  inner  and  outer  eddy  viscosity  formulas 
are  matched  is  also  indicated  (JMATCH) . 

In  Figure  13d,  the  flow  variables  on  the  computational  sublayer  mesh 
adjacent  to  the  lower  boundary  are  given  at  each  station.  The  values  of 
the  normal  distance  (XN)  in  feet  of  each  sublayer  point  are  shown  for  each 
value  of  j  from  1  to  JSL0  =  20,  together  with  the  values  of  the  static 
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temperature  (TEMP)  in  °R,  sublayer  velocity  (U8L)  in  ft/sec,  turbulent  eddy 
viscosity  (EPSLN)  in  lbf-sec/ft‘  and  molecular  viscosity  (VISCOSITY)  in 
lbf-sec/ft2.  For  brevity,  only  those  results  for  i  =  1  to  4  are  shown.  A 
similar  output  follows  (not  shown)  for  the  flow  variables  on  the  computa¬ 
tional  sublayer  mesh  adjacent  to  the  upper  boundary. 

In  Figure  13e,  the  values  of  several  useful  flow  variables  are  shown 
at  each  station  on  the  lower  and  upper  boundaries.  These  variables  are 
p/p  (P/PINF) ,  m  (MASS  BLEED)  in  slug/ft2-sec,  and  c.  (CF).  Stations  ahead 
of  a  leading  edge  have  c^  equal  to  zero.  'Die  values  of  m  are  identical  to 
those  specified  in  the  input  data. 

Finally,  in  Figure  13f  the  values  of  DY  PLUS  and  DYSl,  PLUS  (See  Section 
IV.D.^are  listed  for  the  lower  and  upper  surfaces.  These  values  assist  ■' r. 
determining  whether  the  mesh  distribution  employee  was  satisfactory  in  re¬ 
gards  to  providing  sufficient  flow  field  resolution  (Section  III. A. 2). 

H.  Interpolating  Flow  Variables  at  Mesh  Over  Lap  Using  Program  UPSIRM 

1 .  Introduction 

In  certain  instances  of  mesh  overlapping,  the  heights  of  the  two 
mesh  systems  at  the  restart  station  are  not  equal.  This  is  illustrated  in 
Figure  8,  which  represents  a  typical  system  of  meshes  employed  to  compute 
one  of  the  MCAIR  high  speed  inlet  configurations.  The  flat  plate  portion 
of  the  ramp  ahead  of  the  inlet  entrance  is  computed  using  mesh  A,  whose 
height  is  typically  five  times  the  boundary  Layer  thickness  of  the  restart 
station.  The  height  of  the  overlapping  mesh  B  at  the  restart  station  depend 
on  the  cowl  angle,  and  is  generally  different  from  that  of  mesh  A.  The  progr 
UPSTRM  is  employed  to  interpolate  the  flow  variables  onto  the  upstream  hound 


of  mesh  B  at  the  restart  station.  In  addition,  the  program  provides  initial 
values  for  the  flow  variables  at  all  other  points.  In  this  manner,  a  variety 
of  different  cowl  angles  can  be  considered  for  given  freestream  conditions 
without  the  need  for  recomputing  the  flow  in  Region  A. 

The  program  UPSTRM  is  employed  only  in  the  case  described  above.  In 
those  instances  where  the  heights  of  the  mesh  regions  at  the  restart  station 
are  identical  (e.g.,  meshes  B  and  C  in  Fig.  8),  the  program  UPSTRM  is  not  used. 

The  following  requirements  must  be  satisfied: 

a)  The  restart  station  is  a  vertical  line  in  both  the  upstream  and 
downstream  meshes.  In  its  present  configuration,  the  program  COORD 
insures  that  this  requirement  is  satisfied. 

b)  There  is  a  boundary  layer  on  the  lower  surface  only  in  the  upstream 
region  (mesh  A)  as  illustrated  in  Figure  8. 

c)  The  height  of  the  first  row  of  ordinary  mesh  points  above  the  lower 
surface  at  the  restart  station  is  the  same  in  the  upstream  and  down¬ 
stream  meshes.  This  requirement  must  be  met  in  defining  the  ?  =  0 
boundary  of  the  downstream  mesh  using  program  BNDRY.  All  other 
ordinary  points  maybe  distributed  as  desired  on  this  boundary,  with 
care  taken  to  provide  adequate  resolution  as  discussed  in  Section  III. A. 

The  following  table  provides  information  on  the  resources  required  to 
execute  program  UPSTRM  on  the  CYBER  175. 

TABLE  12.  RESOURCES  REQUIRED  FOR  EXECUTION  OF  PROGRAM  UPSTRM 


Resource 


Details 


Computer  Time 
Input/Output  Time 
Core  Memory 
Files  Required 


<  5  sec  (typical) 

<  5  sec  (typical 

<  150,000  words  (octal) 

INPUT,  OUPUT,  KLOWUP ,  FLOWDN, 
MESHUP,  MESHDN  (See  Section  IV.  II.  5) 
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2.  Description  of  Input  Variables 


The  input  to  program  UPSTRM  consists  of  two  cards.  The  definition 
and  format  of  the  input  data  are  indicated  in  the  required  order  in  fable 
13. 


TABLE  13.  INPUT  DATA  FOR  PROGRAM  UPSTRM 

Line  1:  JSLOLD,  JLOLD,  IL,  JL,  JSL0 ,  JSL1,  JRKF0,  10VLP 

(Format  815) 


For  trail 

Def ini tion 

JSLOLD 

Number  of  points  in  CSL  on 

H  =  0  in  upstream  mosh 

JLOLD 

Number  of  points  in  n-nireet ion 
in  upstream  mesh 

IL 

Number  of  points  in  ^-direction 
in  downstream  mesh 

Ran^.e  or_  Value 
3;  JSLOLliO' 

M  >LO| 


JL  Number  of  points  in  n-diroction  i  <5 

in  downstream  mesh 

JSLI^  Number  of  points  in  CSL  on  3-  ,i:i  vi;- JO 

H  =  0  on  downs  tream  mesh 

JSL1  Number  of  points  in  CSL  on  '  jS  .i-  JO 

H  =  1  on  downstream  mesh 


JREF0 


Value  of  JREF0  in  downstream  nie-ii 


on 


10VL1 


Value  of  i  in  upstream  i'cu'on 
that  defines  restart  station 


1 


TABLE  13.  CONT'D 


Line  2:  UINF,  VINF,  EIINF,  PINF,  EPSINF,  DELTA 


(Format  6F10.A) 


Fortran 

UINF 

VINF 

EIINF 

PINF 

EPSINF 


Def i ni t ion 


U 

oo 

(vertical 
velocity  at  edge 
of  boundary  layer 
at  restart  station) 


Range 


Uni  ts 
Units 


Units 

Uni  ts 


£  at  edge  of  boundary  Units: 
layer  at  restart 
station 


Boundary  layer 
thickness  at 
restart  station 


or  Value 
ft/sec 
f  t/sec 

ft2/sec2 
lbf/ft2 
lbf-sec/f t 


DELTA 


Units:  ft 


Flow  Chart:  Program  UPSTRM 


The  pertinent  features  of  the  program  UPSTRM  are  indicat 


BIT,  IN 


Read  input  data  from 
file  INPUT 


Read  upstream  coordinate 
values  for  y.  .  from  file  MUSUUP 

’  t,j  _ 


Read  downstream  coordinate 

values  for  y.  .  from  file  MF.S'lJiN 

_  i,J . . 


Read  upstream  flow  variables 
from  file  FbOWUP 


Interpolate  variables  in  eomput at inna i 
sublayer  on  lower  surface 


Set  freestream  values  for  variables 
in  computational  sublayer  on  upper  sm far 


Interpolate  variables  in  ordinary  ! 
mesh  at  restart  station 


Provide  initial  values  for  variable 
in  ordinary  mesh  for  r.  0 


Write  flow  data  on  file  FI.OWDN 


Print  interpolated  variables  at  upstream  houndar 


ed  below 


4 .  Source  Code  Notation 


The  definition  of  the  pertinent  source  code  variables  is  given  below. 
The  labeling  conventions  for  the  subscripts  I,  J  and  K  are  the  same  as  in 
Section  IV. D. 


Fortran 

Definition 

Range  or  Value 

DELTA 

Boundary  layer  thickness  at 
restart  station 

Units:  feet 

EI(I,J) 

e.  at  (x.  . ,y .  .) 

1  i,j  i,j 

2  2 

Units:  ft  /sec 

El  INF 

e_-  =  c  T 

v  <» 

2  2 

Units:  ft  /sec 

EPS INF 

c  at  edge  of  boundary  layer 

2 

Units:  lbf-sec/ft 

EPSLN(I, J) 

e .  .  in  OM 
i»  j 

2 

Units:  lbf-sec/ft 

epslsl(i,j,k) 

e.  .  in  CSL 
i»J 

Units:  lbf-sec/ft^ 

EREF(J) 

e  (s  ,n)  in  OM 
eq  o 

Units:  lbf-sec/ft^ 

ESLREF(J.K) 

e  (s  ,n)  in  CSL 
eq  o 

2 

Units:  lbf-sec/ft 

FAC(I,K) 

i/]x  |p  /26p 

1  w 1  w  w 

Units:  ft  ^ 

IL 

Number  of  points  in  4 -direction 
in  downstream  mesh 

<40 

I0VLP 

Value  of  i  in  upstream  region 
that  defines  restart  station 

See  IV. D. 

JL 

Number  of  points  in  n-direction 
in  downstream  mesh 

9  £  JL  £  48 

JLOLD 

Number  of  points  in  n-direction 
in  upstream  mesh 

9  <  JL  £  48 

JSLOLD 

Number  of  points  in  CSL  on 
n  =  0  in  upstream  mesh 

3  <  JSLOLD  <_  20 

JSL0 

Number  of  points  in  CSL  on 
n  =  0  in  downstream  mesh  (need 
not  be  same  as  JSLOLD) 

3  £  JSL0  £  20 

JSL1 

Number  of  points  in  CSL  on 
n  =  1  in  downstream  mesh 

3  £  JSL1  £  20 
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Fortran 

Definition 

Range  or  Value 

JREF0 

Value  of  JREF0  in  downstream 
mesh 

See  IV. D 

P(I,J) 

Pi,i 

2 

Units:  lbf/ft 

PINF 

Poo 

Units:  lbf/ft^ 

QYWL(I,K) 

Heat  transfer  at  wall  in  CSL 

0  (adiabatic  wall) 

QY2(I,K) 

Heat  transfer  normal  to  wall 
at  matching  point 

Units:  lbf/ft-sec 

RHO  ( I ,  J ,  L  ) 

Pi.J 

See  IV.  D 

RHOE(I,J,L) 

Pei,j 

See  IV.  D 

RHOU(I, J,L) 

PUid 

See  IV.  D 

RHOV(I,J,L) 

pv.  . 

See  iv.  I) 

TAUWL  ( I ,  K) 

Wall  shear  stress 

See  IV.  J) 

TAU2(I,K) 

Shear  stress  at  matching  point 

See  IV.  I) 

TSL(I, J,K) 

Static  temperature  in  CSL 

See  IV.D 

0(1, J) 

u.  ,  in  OM 

1  »  J 

S  ee  l  V .  I) 

USL(I, J,K) 

u  in  CSL 

See  IV. I) 

UINF 

U 

oo 

Units:  ft/sec 

V(I,J) 

v.  .  in  OM 

1,  j 

Units:  ft/sec 

VINF 

V 

GO 

Units:  ft/sec 

XMUSL(I, J,K) 

Vi  in  CSL 

See  CV.I) 

YUP(I, J) 

y.  .  in  upstream  mesh 
i » J 

Units:  feet 

YDOWN(I, J) 

y.  ,  in  downstream  mesh 

i>  j 

U  n i L  s :  feet 

YSLUP(J) 


Mesh  point  distribution  in  CSL 
at  restart  station 


f  ee  t 


5. 


File  Structure 


The  program  employs  four  files,  in  addition  to  the  conventional  files 
INPUT  and  OUTPUT.  The  descriptions  are  indicated  below.  For  further 
information  see  Section  IV. E. 


File 

Description 

Data  Structure 

FLOWUP 

Input  file  of  flow  variables  in 
upstream  region 

See  file  RSTART 

FLOWDN 

Output  file  of  flow  variables 
in  downstream  region 

Same  as  FLOWUP 

MESHUP 

Input  file  of  coordinate  trans¬ 
formation  data  for  upstream 
region 

See  file  XY  in 
Section  IV. D 

MESHDN 

Input  file  of  coordinate  trans¬ 
formation  data  for  downstream 
region 

Same  as  FLOWDN 

All  four  files  are  required  for  each  execution  of  the  program. 

6.  Output 

The  printed  output  consists  of  (1)  the  interpolated  flow  variables 
u,  v,  e^  p  and  e  on  the  ordinary  mesh  of  the  downstream  region  at  the 
restart  station,  and  (2)  the  interpolated  flow  variables  T,  u,  £  and  y  on 
the  computational  sublayer  mesh  (r)  =  0  boundary)  of  the  downstream  region 
at  the  restart  station.  The  format  is  similar  to  that  employed  by  the 
Navier-Stokes  code  (see  Figures  13b  and  13d) .  It  should  be  noted  that  the 
messages  "RM  ERROR  0142  on  LFN  MESHUP"  and  "RM  ERROR  0142  ON  LFN  MESHDN" 
may  appear  in  the  dayfile.  These  messages  are  of  no  significance,  and  may 
be  ignored. 
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Region  4 


Region  3 

Region  2 

Region  1 


Figure  4.  Mesh  Splitting 


•  Computational 
Sublayer  Mesh 


a)  Computational  Sublayer  on  1  =  0 


b)  Computational  Sublayer  on  p  =  1 


Figure  7.  Computational  Sublayer  Geometry 
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Figure  10.  Sample  Output:  Program  1SNURY 
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Figure  11a.  Coordinate  Transformation  for  MCAIR  inlet  (Upstream) — Case  35 
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Figure  12a.  Sample  Output:  Program  COORD 
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Figure  12e.  Sample  Output:  Program  COORD 
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